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1.0  SUMMARY 


A  one-dimensional  software  tool  designed  for  laboratory  modeling  (NUMIT)  has  been  refined, 
developed,  and  greatly  enhanced  so  that  it  can  be  used  to  model  and  hopefully  predict  the  deep- 
dielectric  charging  of  dielectrics  on  spacecraft.  The  new  tool,  called  AF-NUMIT2,  can  now 
model  isotropically  incident  electron  fluxes  from  realistic  energy  spectra  that  change  with  time. 
The  model  has  been  extended  to  incident  electron  fluxes  an  order  of  magnitude  lower  than  before 
(from  >100  keV  to  >  10  keV),  and  aspects  of  the  model  have  been  checked  against  Monte  Carlo 
simulations.  In  order  to  verify  the  utility  of  AF-NUMIT2,  dozens  of  simulated  runs  have  been 
made  using  realistic,  isotropic  electrons  incident  on  Kapton,  representing  a  wide  range  of  actual 
energy  spectrum  fluxes  measured  by  CRRES.  Kapton  has  been  modeled  with  both  dark 
conductivities  and  radiation  induced  conductivities  extending  over  two  orders-of-magnitude.  The 
result  has  been  a  surprisingly  consistent  relationship  between  maximum  internal  electric  field 
strengths  and  conductivity. 


2.0  INTRODUCTION 

The  original  NUMIT,  short  for  numerical  iteration,  was  designed  and  described  by  A.  R. 
Frederickson.  [1]  [2]  It  was  a  one-dimensional  model  intended  primarily  to  aid  in  the 
interpretation  of  laboratory  testing  during  which  mono-energetic  electron  beams  had  normal 
incidence  on  various  dielectric  slabs.  The  basic  concepts  of  NUMIT  have  been  utilized  in  the 
present  work  and  must  be  understood  before  the  development  of  AF-NUMIT2  can  be  explained. 
NUMIT  necessarily  combines  three  different  algorithms  or  models:  energy  deposition,  electron 
deposition,  and  charge  transport.  Summarized  here  are  the  primary  algorithms  used  when  the 
author  worked  with  Frederickson  in  2000: 

1 .  EDEPOS  An  algorithm  is  required  to  determine  the  deposition  of  energy  as  a  function  of 
depth  by  the  incident,  high-energy  electrons.  NUMIT  accomplishes  this  by  utilizing  the 
well-known  EDEPOS  algorithm,  [3]  which  in  fewer  than  fifty  lines  of  code  compiles  the 
results  of  numerous  experiments  and  Monte  Carlo  simulations  for  various  materials  and 
incident  electron  energies. 

2.  FredBell  An  algorithm  is  required  to  determine  where  the  incident  electrons  stop  in  the 
dielectric  as  a  function  of  depth.  We  use  an  algorithm  by  Frederickson,  Bell,  and  Beidl 
[4]  that  draws  heavily  on  EDEPOS  and  incorporates  Monte  Carlo  simulations.  For 
convenience,  we  name  this  algorithm  “FredBell.” 

3.  Transport  Model  A  transport  model  is  required  to  determine  how  the  charge  moves 
about  within  the  dielectric.  NUMIT  detennines  the  conductivity  in  the  dielectric  that 
results  from  the  energy  deposition  by  utilizing  the  well-known  Fowler  Model  [5]  for 
radiation  induced  conductivity  (RIC).  Then,  the  current  can  be  predicted  once 
electrostatics  is  used  to  calculate  the  electric  field  resulting  from  the  charge  deposition. 
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Recent  work  has  shown  that  internal  electric  fields  are  the  best  predictor  of  deep-dielectric 
pulsing.  [6]  [7]  Therefore,  it  is  the  calculation  of  the  electric  field  internal  to  the  dielectric,  and 
the  prediction  of  how  it  might  change  in  time,  that  makes  NUMIT  a  potentially  useful  approach 
for  assisting  the  spacecraft  designer. 

Frederickson’s  original  algorithm  has  been  greatly  expanded  and  enhanced  during  the  course  of 
the  present  contract  in  order  to  enable  NUMIT  to  model  what  might  be  happening  to  dielectrics 
in  space.  The  development  of  a  limited,  laboratory  modeling  code  into  a  more  sophisticated  code 
AF-NUMIT2  that  is  useful  for  predicting  the  impact  of  the  space  environment  on  Air  Force 
systems  is  a  multifaceted  task.  The  goals,  as  delineated  in  the  SOW,  are  listed  here  with  only  the 
order  changed  to  reflect  the  sequence  in  which  they  are  addressed  in  the  present  report: 

1.  (Sec  3.2)  Establishment  of  a  definitive,  controlled  source  code  for  AF-NUMIT2 

2.  (Sec  3.3)  Method  for  accounting  for  electrons  incident  at  different  angles 

3.  (Sec  3.4)  Provision  for  arbitrary  electron  energy  spectrums 

4.  (Sec  3.5)  Provisions  for  time-sequenced  incident  electron  flux  environments 

5.  (Sec  3.6)  Extension  of  model  for  incident  electrons  to  lower  energy  range 

6.  (Sec  3.7)  Verification  of  the  energy  deposition  model  in  NUMIT 

7.  (Sec  3.8)  Development  of  time-delayed  and  temperature-dependent  radiation  induced 
conductivity  algorithms 

A  detailed  discussion  of  each  of  these  goals  will  be  in  the  next  section. 


3.0  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 

3.1  The  Original  NUMIT  Algorithms 

The  plan  for  the  work  reported  here  was  to  lean  heavily  on  the  original  methods  and  assumptions 
utilized  by  Frederickson  in  his  NUMIT  simulation  in  order  to  devise  a  method  of  predicting  the 
electric  fields  that  develop  within  spacecraft  dielectrics  subjected  to  the  space  environment. 
Frederickson  utilized  three  algorithms  for  his  one-dimensional  laboratory  model,  each  of  which 
have  been  modified  as  part  of  the  development  of  AF-NUMIT2.  Accordingly,  each  will  be 
described  briefly  in  this  section. 

3.1.1  Energy  Deposition  Profile  (EDEPOS).  EDEPOS  (short  for  energy  deposition)  was 
created  in  1974  by  Tabata  and  Ito  [3]  and  has  become  an  important  and  widely  used  algorithm 
for  determining  the  profile  of  energy  deposition  by  normally  incident  electrons.  EDEPOS  uses 
the  results  of  more  than  twenty  experimental  studies  and  more  than  a  half-dozen  Monte  Carlo 
method  calculations.  At  the  time,  the  leading  algorithm  had  been  produced  by  Kobetich  and  Katz 
[8],  but  Tabata  and  Ito  demonstrated  that  their  version  was  valid  over  a  wider  range  of  incident 
energies  with  greater  accuracy. 

In  addition  to  the  energy  of  the  incident  mono-energetic  electrons,  EDEPOS  requires  only  the 
input  of  the  insulator's  atomic  number,  atomic  weight,  and  density.  Effective  atomic  numbers  can 
be  used  so  that  the  range  of  applicability  is  from  about  5.3  to  82.  The  incident  electron  energy 
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range  is  0.1-20  MeV.  The  algorithm  was  published  as  only  forty-five  lines  of  FORTRAN  code. 
[3]  Its  ready  availability,  rapid  computation  time,  wide  applicability,  and  ease  of  use  has  made 
EDEPOS  a  tool  of  choice  for  many  modelers.  The  output  of  EDEPOS,  when  multiplied  by  the 
incident  electron  beam  current  Jmc,  is  the  energy  deposition  profiles,  the  dose  rate  D(x),  required 
to  determine  RIC  in  NUMIT. 

3.1.2  Incident  Electron  Deposition  Profile  (FredBell).  In  order  to  model  the  movement  of 
charge  within  the  dielectric,  one  must  know  where  the  high-energy  incident  electrons  stop.  In 
1995,  Frederickson,  Bell,  and  Beidl  [4]  produced  an  algorithm  that  superseded  all  earlier  guesses 
at  where  the  charge  was  being  deposited.  Frederickson  et  al.  took  Tabata's  EDEPOS  and 
modified  it,  using  equations  that  “crudely”  simulated  physical  processes.  These  equations  were 
then  assembled  into  an  algorithm  with  six  fitting  parameters.  By  comparison  with  Monte  Carlo 
simulations,  Frederickson  was  able  to  determine  analytical  functions  for  five  of  the  six 
parameters.  He  left  the  sixth  (called  “5”)  to  be  detennined  from  a  six-column  table. 

Because  it  draws  on  the  energy  deposition  algorithm  by  Tabata,  the  FredBell  algorithm  requires 
use  of  EDEPOS — it  cannot  be  used  as  a  stand-alone  algorithm.  The  required  inputs  are  the  same 
as  for  EDEPOS,  and  the  output  is  the  fraction  of  incident  electron  flux  as  a  function  of  depth. 
When  multiplied  by  the  incident  electron  beam  current  Jmc,  FredBell  provides  the  initial  electron 
current  profile  Jo(x)  that  is  required  in  NUMIT. 

3.1.3  Charge  Transport  Algorithm.  Once  the  energy  and  electron  deposition  are  detennined  as 
a  function  of  depth  in  the  dielectric,  NUMIT  creates  a  transport  algorithm  based  on  electrostatics 
and  Fowler's  conductivity  model.  [5]  Although  NUMIT  can  be  modified  for  different 
configurations,  the  simplest  setup  is  illustrated  in  Figure  1.  A  dielectric  slab  is  sandwiched 
between  two  electrodes  that  can  be  held  at  the  same  electric  potential  or  at  a  fixed  potential 
difference. 

In  Figure  1,  we  take  the  left  side  of  the  dielectric  to  be  x  —  0  and  measure  the  depth  from  that 
point.  The  incident  high-energy  electrons  are  assumed  to  come  from  the  left  and  penetrate  into 
the  dielectric.  For  the  purposes  of  modeling,  the  left  electrode  can  be  taken  to  be  immeasurably 
thin,  such  as  conductive  paint.  Initially,  in  the  absence  of  an  accumulation  of  charge,  the  incident 
electrons  create  a  position-dependent  current  density  J0(x). 

The  initial  current  density  is  virtually  always  a  monotonically  decreasing  function  of  depth  x 
because  the  electrons  are  slowed  and  absorbed  by  the  dielectric,  and  significantly  fewer  electrons 
reach  the  far  side  than  enter  the  dielectric.  As  the  electrons  are  absorbed,  they  give  up  their 
energy  to  the  dielectric.  The  rate  of  this  energy  deposition  is  called  the  dose  rate  D(x)  and  is 
given  in  rad/sec.  Unlike  Jo,  the  dose  rate  is  not  a  monotonically  decreasing  function  of  x,  but 
typically  peaks  before  exhibiting  a  tail.  In  NUMIT,  D(x)  and  J0(x)  are  approximated  by  the 
EDEPOS  and  FredBell  algorithms  discussed  above. 
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Figure  1.  Schematic  Diagram  of  Situation  Modeled  by  NUMIT 


As  shown  in  Figure  1,  the  left  electrode  at  x  =  0  is  grounded,  while  the  right  electrode  at  x  —  L 
is  held  at  a  constant  potential  V,  which  very  often  is  set  to  zero.  The  dielectric  is  assumed  to  have 
a  constant  permittivity  e.  Its  thickness  L  is  assumed  to  be  small  compared  to  its  other  dimensions 
so  that  the  problem  can  be  treated  in  one  dimension.  We  wish  to  find  the  charge  density  p,  the 
electric  field  E,  and  the  current  density  J  as  functions  of  position  x  and  time  t. 

There  are  three  major  equations  required  to  relate  p,  E,  and  J.  The  first  two  are  the  continuity 
equation  and  the  differential  fonn  of  Gauss's  Law.  In  one-dimension,  these  are 

dj(x,  t)  dp(x,t) 

dx  dt 


dE(x,t )  p(x,t ) 
dx  e 

(2) 


The  third  equation  required  represents  the  transport  of  charge.  Utilizing  the  Fowler  Model  for 
radiation  induced  conductivity,  we  see  that  for  the  same  electric  field,  the  current  will  be  higher 
when  the  dose  rate  is  greater: 

J(x,t)  =/0W  +  [g0  +  kb{x)]E{x,t) 

(3) 

Here  g 0  is  the  dark  conductivity,  k  is  the  coefficient  of  radiation-induced  conductivity  (RIC),  and 
D(x)  is  the  dose  rate.  The  dark  conductivity  is  the  conductivity  of  the  insulator  when  it  is  not 
radiated.  Because  it  is  an  insulator,  the  dark  conductivity  is  extremely  low  and  difficult  to 
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determine  accurately  for  many  important  materials — particularly  in  space  applications.  [9]  [10] 

[11]  [12] 


Equations  (l)-(3)  can  be  used  to  describe  the  transport  of  charge  in  the  dielectric  by  taking 
successive  time  steps.  The  initial  current  Jo(x)  is  due  to  the  incident  high-energy  electrons  and  is 
given  by  the  FredBell  algorithm.  Then  the  Continuity  Equation  (1)  enables  us  to  find  the  charge 
deposited  p(x.  At)  a  moment  At  in  time  later.  Once  we  know  p(x,  At),  Gauss's  Law  (2)  enables  us 
to  find  E(x,  At).  Combining  our  knowledge  of  E(x,  At)  with  the  dose  rate  D(x)  given  by  the 
EDEPOS  algorithm.  Equation  (3)  enables  us  to  detennine  the  movement  of  the  deposited  charge. 
Of  course,  electrons  continue  to  be  incident  on  the  material  so  Jo(x)  must  be  added  to  find  the 
new  total  current  J(x,  At)  within  the  dielectric.  This  cyclical  time  stepping  process  can  continue 
as  long  as  desired  or  until  equilibrium  is  achieved. 

The  numerical  algorithm  requires  careful  choice  of  the  sizes  of  the  spatial  increment  Ar  across 
the  dielectric  and  the  time  step  size  At.  The  choice  of  Ax  detennines  the  resolution  of  all  the 
calculated  values,  but  too  small  a  value  will  make  the  simulation  take  longer  than  necessary. 

Each  Ax  can  be  thought  of  as  a  bin  which  holds  the  charge  in  that  region  of  the  dielectric.  The 
charge  is  considered  to  reside  in  the  center  of  each  bin,  called  a  node.  For  normally  incident 
electrons,  the  situation  for  which  NUMIT  was  originally  designed,  30-200  nodes  usually  provide 
adequate  resolution.  However,  when  isotropically  incident  electrons  are  included,  the  number  of 
nodes  should  normally  be  increased  because  deposition  will  trend  shallower.  In  AF-NUMIT2, 
even  when  considering  electrons  incident  at  angles  other  than  perpendicular  to  the  surface,  the 
spatial  nodes  continue  to  represent  the  perpendicular  depth  from  the  surface. 

Similarly,  the  time  step  At  must  be  chosen  very  carefully.  If  it  is  too  small,  the  program  runs  for 
an  excessively  long  time,  diminishing  its  usefulness.  If  it  is  too  large,  the  net  charge  begins  to 
fluctuate  wildly.  The  fluctuations  occur  because  if  a  time  step  is  too  long,  then  the  algorithm  will 
calculate  more  charge  having  been  removed  from  the  node  than  the  corresponding  bin  actually 
holds. 

Additional  details  on  NUMIT,  including  updated  calculations  of  surface  charging  on  the 
electrodes,  can  be  found  elsewhere. 1  [13] 


3.2  Establishment  of  a  Definitive,  Controlled  Source  Code  for  AF-NUMIT2 

The  original  NUMIT  was  written  by  Frederickson  in  FORTRAN.  Input  and  output  was  very 
tedious  and  required  the  use  of  another  software  application  for  decent  plots  of  the  results.  When 
the  author  worked  with  Frederickson  in  2000,  and  inserted  the  FredBell  algorithm  for  the  first 


1  See  also  Brian  P.  Beecken,  “Development  of  the  NUMIT  Simulation  for  Modeling  Deep-dielectric  Charging  in  the 
Space  Environment,"  a  report  to  AFRL/RVBX,  August  13,  2009,  particularly  Appendix  A. 
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time,  the  FORTRAN  code  had  about  half  its  lines  and  subroutines  commented  out.  In  2007,  the 
group  at  JPL  responsible  for  continuing  Frederickson’s  work  after  his  death  wrote  “NUMIT  has 
not  been  configuration  controlled.”  [14] 

At  the  time  of  his  inclusion  of  the  FredBell  algorithm,  the  author  discovered  that  the  way  it 
handled  delta  ray  current"  could  not  be  utilized  without  double  precision  in  FORTRAN.  For 
expediency,  in  consultation  with  Frederickson,  it  was  decided  to  leave  out  delta  rays  because 
they  were  a  relatively  small  effect.  FredBell  was  also  limited  in  usefulness  by  the  requirement  of 
looking  up  a  parameter  called  “B”  in  a  published  table.  Each  change  of  material  and/or  energy 
required  finding  the  appropriate  value  for  B  in  the  table  and  inserting  it  manually  in  the  code. 
These  limitations  have  recently  been  independently  documented  and  corrected  at  JPL.  [15] 

AF-NUMIT2  is  written  in  MATLAB™  which  allows  for  easy  input  of  relevant  parameters  and 
the  rapid  output  of  publication  quality  graphs,  and  it  is  inherently  comparable  to  double  precision 
in  FORTRAN.  The  subroutines  EDEPOS  and  FredBell  have  been  adjusted  and  combined  into 
one  MATLAB  function  named  BeeckenTabFred.  The  code  accounts  for  delta  ray  current  and  has 
been  revised  for  efficiency  so  that  it  runs  significantly  faster.  Speed  improvements  are  important 
because  AF-NUMIT2  requires  many  loops  through  this  function  as  it  simulates  the  space 
environment  with  isotropically  incident  electrons.  A  short  function  called  “getparameterB”  was 
written  to  interpolate  the  Frederickson  table  for  “B”,  and  it  allows  the  program  to  run  without 
additional  user  inputs. 

Figure  2  is  a  flow  chart  for  AF-NUMIT2.  The  blue  Time  Loop  portion  is  essentially  the  original 
NUMIT  as  described  above,  written  in  MATLAB  with  the  aforementioned  enhancements.  The 
yellow  Spectrum  Loop  block  represents  the  bulk  of  the  additions  that  had  to  be  made  to  the  code 
in  order  to  predict  charging  in  realistic  space  environments.  The  green  Inputs  block  describes  the 
far  more  complicated  inputs  that  are  required  for  the  program  to  model  the  space  environment, 
and  the  red  Plotting  block  illustrates  the  new  standard  outputs. 


2  Delta  ray  current  is  a  name  for  electrons  that  are  knocked  loose  in  the  material  by  the  high-energy  incident 
electrons.  On  average  they  move  deeper  within  the  material. 
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Figure  2.  Flow  Chart  for  AF-NUMIT2 

AF-NUMIT2  has  been  produced  in  two  versions.  The  first  is  written  as  a  standard  Matlab™ 
program  with  extensive  commenting.  It  is  purposely  designed  to  be  easily  modified  and 
improved  by  users  who  have  some  familiarity  with  Matlab.  Instructions  on  its  use,  the  first  100 
lines  of  code  necessary  for  inputs,  and  a  list  explaining  the  variables  used  within  the  code  (to  aid 
in  modification),  have  been  provided  in  Appendices  A,  B,  and  C.  Calculation  of  effective  atomic 
number  and  weight  for  composite  materials  is  explained  in  Appendix  D.  The  second  version  of 
AF-NUMIT2  has  a  Graphical  User  Interface.  It  contains  the  same  algorithms,  but  should  be  an 
easier  version  of  AF-NUMIT2  to  use  initially.  Unfortunately,  it  runs  much  slower  and  currently 
is  somewhat  flaky.  Therefore,  as  of  this  writing,  that  version  is  not  being  submitted,  nor  was  it 
included  in  the  SOW.  There  has  been  no  provision  to  copyright  AF-NUMIT2  or  otherwise  give  it 
configuration  control. 
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3.3  Method  for  Accounting  for  Electrons  Incident  at  Different  Angles 


A  primary  limitation  in  the  conversion  of  NUMIT  from  a  laboratory  tool  to  a  model  useful  for 
spacecraft  design  is  not  that  NUMIT  is  a  one-dimensional  model  of  charge  transport,  but  rather  it 
is  NUMIT ’s  inherent  use  of  an  electron  beam  that  strikes  the  dielectric  surface  at  only  normal 
incidence.  Many  dielectrics  of  interest  on  spacecraft  are  essentially  flat,  so  the  one-dimensional 
model  of  charge  transport  within  the  dielectric  is  an  appropriate  approximation.  The  issue  then  is 
what  is  the  appropriate  way  to  handle  an  isotropic  flux  of  electrons  striking  the  surface? 

3.3.1  Determination  of  Electron  Flux  as  a  Function  of  Incident  Angle.  The  determination  of 
the  correct  electron  flux  to  use  is  non-trivial.  The  space  environment  nominally  has  an  isotropic 
electron  flux.  For  a  particular  angle  from  the  normal  6,  we  need  to  find  the  appropriate  fraction 
Jmc(9)  of  the  total  flux  that  will  strike  the  dielectric  surface  and  then  penetrate.  What  is  the 
appropriate  incident  flux  as  a  function  of  angle? 

Getting  the  correct  angle  dependent  incident  electron  flux  is  absolutely  crucial  to  a  realistic 
prediction  of  spacecraft  charging.  Accordingly,  much  effort  has  been  expended,  and  two 
approaches  have  been  considered  and  discarded.  At  first,  it  might  seem  reasonable  to  assume  that 
the  calculation  is  merely  the  usual  determination  of  the  number  of  electrons  that  strike  the  top 
surface  of  a  disk  in  an  isotropic  flux.  Such  a  calculation  is  well-known,  but  not  relevant  because 
we  are  interested  in  penetrating  flux.  Secondly,  it  might  seem  that  we  need  to  know  the  amount 
of  charge  that  penetrates  and  then  passes  through  the  imaginary  line  drawn  normal  to  the  surface 
along  which  the  one-dimensional  charge  transport  algorithm  will  proceed.  Significant  work  was 
done  on  this  possible  approach  to  calculate  the  appropriate  Jinc(0)  for  such  a  situation,  and  it  was 
reported  in  two  published  proceedings  papers.  [13]  [16]  Unfortunately,  this  approach  also 
appears  now  to  be  incorrect.  The  simplest  explanation  for  the  error  is  conceptual.  Because  an 
imaginary  line  normal  to  the  surface  has  no  thickness,  a  calculation  of  particle  flux  passing 
through  it  can  have  no  physical  significance. 

The  correct  approach  is  to  consider  the  flat  dielectric  as  being  large  enough  that  there  are  no  edge 
effects.  With  the  incident  electron  flux  being  uniformly  isotropic,  the  deposition  of  charge  within 
the  material  will  be  a  function  of  depth  uniformly  throughout  the  material.  In  other  words,  for  a 
given  depth,  the  layer  of  deposited  charge  at  that  depth  is  uniform. 

Going  back  to  the  calculation  of  the  number  of  electrons  incident  on  the  surface  of  a  disk  to 
determine  how  the  incident  flux  depends  on  angle,  however,  is  incorrect  because  it  treats  the 
electrons  as  electromagnetic  radiation.  A  radiation  calculation  thinks  in  terms  of  the 
effective  area  of  the  detector  surface  as  a  function  of  angle — one  that  shrinks  as  the  angle  of 
incidence,  measured  from  the  normal,  increases.  However,  the  implicit,  but  usually  unstated, 
assumption  is  non-penetration.  For  penetrating  electrons,  see  Figure  3,  there  would  be  no 
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such  thing  as  effective  area  as  all  electrons  that  are  incident  go 
into  the  material  regardless  of  angle  (neglecting  backscatter 
effects).  This  change  in  treatment  required  much  thought,  many 
trial  calculations,  and  was  considered  at  great  length  with 
members  of  AFRL/RVBX.  The  result  is  stunningly  simple:  for 
isotropic  flux,  the  incident  flux  is  the  same  at  all  angles.  For 
AF-NUMIT2,  the  calculation  of  the  electron  current  per  unit 
area  Jinc  is  done  simply  by  taking  the  incident  flux  /incident  given  in 


Figure  3.  Illustration  of 
Penetrating  Model 


units  of  ( 1  /  cm  sr  sec)  at  a  specific  energy  and  multiplying  by  the 
steradians  in  one  hemisphere  so  as  to  include  all  electrons  incident  on  the  material: 


J  inc  27r/incident. 

If  this  treatment  is  correct,  the  incident  flux  does  not  depend  on  angle.4  However,  for  the 
purposes  of  AF-NUMIT2,  electrons  striking  at  different  angles  must  be  treated  differently 
because  they  will  have  different  penetration  depths.  Therefore,  the  model  must  keep  track  of 
incident  angles  and  this  is  done  by  dividing  the  flux  into  angular  bins  of  size  2vp  (measured  in 
degrees)  at  each  angle  0.  The  result  is 

2  ip 

J inc(^)  —  90°  ^*nc 

Note  that  the  fraction  must  use  angles  not  steradians  because  the  horizontal  direction  from  which 
the  electrons  come  does  not  matter.  We  only  want  the  fraction  of  all  incident  electrons  at  a 
certain  angle  0  to  the  normal. 

3.3.2  Determination  of  Energy  and  Electron  Deposition  Profiles  as  Functions  of  Angle.  We 

make  the  assumption  that  most  dielectrics  of  interest  are  going  to  be  isotropically  homogeneous 
(certainly  not  crystalline).  Of  course  many  dielectrics  will  be  polymers,  but  we  have  to  start 
somewhere.  We  hope  that  the  polymer  chains  are  not  consistently  oriented  in  the  same  direction 
so  that  incident  electrons  will  not  penetrate  significantly  more  in  one  direction  than  another.  If 
homogeneity  is  a  reasonable  approximation,  then  as  the  incident  high-energy  electrons  travel 
through  the  dielectric  in  the  direction  of  their  straight  line  trajectory  (at  a  particular  angle  of 
incidence)  the  depth  of  penetration  and  the  rate  of  energy  deposition  should  be  relatively 


3  We  are  grateful  to  William  (Bob)  Johnston  for  significant  and  helpful  conversations  during  August  of  201 1  and  to 
David  Cooke,  Dale  Ferguson,  and  Adrian  Wheelock  for  a  long  discussion  and  eventual  consensus  on  this  approach 
at  a  meeting  on  25  July  2012. 

4  In  the  AF-NUMIT2  code,  ./incident  is  for  specific  energy  channels,  just  as  it  is  here.  However,  in  the  code,  instead  of 
Jinc  in  terms  of  (1/cm2)  we  actually  use  Jjncid  which  is  (A/cm2). 
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independent  of  the  angle.  We  are  talking  strictly  about  the  situation  after  the  electrons  have 
entered  the  material.  As  we  have  already  shown,  the  number  of  electrons  striking  the  material  is 
believed  to  be  angle  independent.  The  number  of  electrons  that  are  potentially  backscattered, 
however,  is  obviously  angle-dependent  and  is  considered  in  the  next  subsection. 

Now  we  can  use  the  EDEPOS  and  FredBell  algorithms  (whether  separately  or  in  their  combined 
form  BeeckenTabFred)  to  approximate  energy  and  charge  depositions  for  any  particular  angle  of 
incidence  9.  On  average,  electrons  will  travel  the  same  distance  within  the  material,  however  that 
distance  must  be  scaled  back  to  a  perpendicular  depth — a  simple  trigonometric  correction.  If  we 
call  x  the  perpendicular  distance  into  the  material,  what  NUMIT  treats  as  depth,  then  the 
distances  along  the  path  at  an  angle  of  incidence  0  given  by  BeeckenTabFred  can  be  called  x@. 
The  appropriate  conversion  becomes: 


Xg  COS  6  —  X 

The  BeeckenTabFred  algorithm  will  return  to  AF-NUMIT2,  at  various  perpendicular  depths 
measured  in  x,  the  values  for  the  energy  deposited  per  electron  (EDEPOS)  and  the  fraction  of 
electrons  (Jofrac)  that  continue  on  in  the  chosen  direction  detennined  by  the  incident  angle.  Both 
values  must  be  multiplied  by  a  factor  representing  the  current  density  of  electrons  that  entered 
the  material  at  a  particular  angle  9  from  the  normal.  This  approach  is  clearly  an  approximation. 
However,  if  isotropic  incidence  is  going  to  be  included  in  AF-NUMIT2,  in  the  absence  of  a 
better  approach  we  must  do  something.  The  validity  of  this  approach  has  been  examined  using 
Monte  Carlo  simulations;  it  will  be  discussed  in  Subsection  3.7.4. 

3.3.3  The  Issue  of  Backscatter.  When  electrons  are  incident  on  a  material,  it  is  known  that  some 
of  the  electrons  penetrate  while  others  are  reflected  or  backscattered.  Very  little,  if  any,  data  is 
known  to  exist  for  backscattering  from  dielectric  materials — particularly  as  a  function  of  angle. 
The  problem  is  complicated  by  secondary  electron  yields,  auger  electrons,  and  differing  material 
structures.5  Nevertheless,  there  is  a  real  effect.  Certainly  it  will  be  angle  dependent,  and 
something  must  be  done  to  at  least  make  an  approximation  for  the  sake  of  AF-NUMIT2  charging 
predictions.  The  approach  we  took  is  explained  here. 

Because  of  backscatter,  it  is  necessary  to  correct  the  incident  flux  by  a  transmission  factor  1  -B, 
where  B  is  the  fraction  of  electrons  that  backscatter.  When  Frederickson  developed  the  FredBell 
algorithm  for  electron  current  and  deposition,  [4]  working  at  that  time  strictly  with  nonnal 
incidence,  his  B  term  depended  on  electron  energy  and  the  material’s  atomic  number.  Later, 


5  During  a  meeting  on  15  August  2011  at  AFRL/RVBX,  several  relevant  comments  were  made  such  as: 
“Backscatter  yields  on  dielectrics  might  not  exist,  even  for  normal  incidence.”  “The  Nascap  model  is  not  help  fid.” 
“How  can  such  basic  data  not  be  available?” 
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Frederickson  developed  another  equation  that  included  dependence  on  incident  angle  0,  in 
addition  to  the  incident  electron  energy  T0  and  the  atomic  number  Z: 


B  =  (0.14  -  0.039ro)(0.1Z°-9)  +  Q12 

This  equation  for  backscatter  was  reported  in  Appendix  A.  1.2  in  a  paper  by  Davis.  [17]  It  would 
seem  obvious  that  the  fraction  of  electrons  backscattering  will  change  as  the  angle  of  incidence 
changes,  and  this  equation  seems  to  provide  a  realistic  estimate  of  that  dependence  on  angle. 


The  transmission  factor  should  be  unity  when  the  backscatter  fraction  B  is  zero.  Unfortunately, 
for  close  to  normal  incidence,  this  equation  becomes  negative  for  higher  electron  energies.  In 
such  cases,  using  it  will  result  in  a  transmission  factor  greater  than  unity.  Having  more  electrons 
in  the  simulation  than  are  actually  incident  on  the  dielectric  would  seem  to  reduce  realism.  It  is 
the  first  term  in  B  that  can  go  negative.  Therefore,  we  have  modified  the  equation  in  AF- 
NUMIT2  by  normalizing  the  entire  expression,  but  only  when  the  first  term  goes  negative. 
Typical  results  are  presented  in  Figure  4. 

Close  examination  of  Figure  4  shows  that  at  the  5  MeV  incident  electron  energy  more 
backscatter  occurs  at  lower  atomic  numbers,  whereas  at  1  MeV  less  backscatter  occurs  at  lower 
atomic  numbers.  This  difference  is  inherent  in  the  equation  for  B  and  is  not  caused  by  our 
normalization  at  the  higher  energy.  The  crossover  occurs  at  about  3  MeV.  We  do  not  know  if  this 
is  a  realistic  representation  of  the  actual  situation. 


3.3.4  Putting  it  all  together:  The  Incident  Angle  Loop.  We  have  shown  that  there  are  at  least 
three  steps  necessary  to  handle  isotropic  incidence:  determining  the  electron  flux  as  a  function  of 


Figure  4:  Sample  Plots  of  Transmission  after  Backscatter 


angle,  calculating  the  current  and  energy  deposition  profiles  in  the  dielectric  that  result  from 
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electrons  at  that  angle,  and  allowing  for  a  dependence  on  incident  angle  of  the  backscattered 
electrons.  But  how  is  this  accomplished  in  the  code? 

The  code  must  consider  each  angle  of  incidence  separately.  The  angles  of  incidence  9  are  divided 
into  angular  bins,  the  size  of  these  correspond  to  the  resolution  chosen  by  the  user.  A  typical 
number  of  bins  might  be  18  between  normal  incidence  and  90°,  correlating  to  5°  for  each  angular 
bin  (called  VV)  with  9  centered  within  the  bin.  As  illustrated  in  the  Flow  Chart  of  Figure  2,  a 
loop  is  instituted  for  each  angular  bin.  Each  time  through  this  loop  the  function  BeeckenTabFred 
calculates  the  energy  and  electron  deposition  profiles.  After  these  profiles  are  detennined  for 
each  angle  9,  the  main  AF-NUMIT2  program  adds  the  charge  deposited  in  each  spatial  bin  and 
separately  adds  the  energy  deposited  in  each  spatial  bin.  As  in  the  original  NUMIT,  the  spatial 
bins  are  simply  defined  as  distance  (or  depth)  increments  of  A x  from  the  surface  as  measured 
perpendicularly. 

AF-NUMIT2  will  by  default  use  the  backscatter  function.  Although  its  validity  and  accuracy  is 
uncertain,  it  is  important  to  account  for  backscatter  when  considering  isotropic  scattering. 
Otherwise,  due  to  the  large  number  of  electrons  incident  at  very  shallow  angles  (large  9),  surface 
charging  will  unrealistically  dominate  a  model  designed  to  predict  deep-dielectric  charging. 

3.3.5  Electron  Beams.  AF-NUMIT2  still  does  electron  beams,  either  mono-energetic  or  multi- 
energetic.  The  number  of  angular  bins  must  simply  be  set  to  one.  Other  input  differences 
required  for  beams  are  described  in  Section  3.4. 


3.4  Provision  for  Arbitrary  Electron  Energy  Spectrums 

The  original  NUMIT,  being  designed  for  a  mono-energetic  electron  beam,  did  not  have  provision 
for  the  wide  spectrum  of  energies  that  electrons  incident  in  a  space  environment  carry.  Two 
aspects  of  AF-NUMIT2  had  to  be  developed  to  address  this  situation:  simultaneous  modeling  of 
multiple  electron  energies  and  the  input  of  those  energies. 

3.4.1  Modeling  multiple  electron  energies.  AF-NUMIT2  has  an  energy  loop  as  shown  in  the 
yellow  Spectrum  Loop  block  of  Figure  2.  For  each  designated  incident  electron  energy,  the  code 
utilizes  the  function  BeeckenTabFred  to  calculate  a  depth  profile  of  the  energy  deposited  per 
electron  and  the  fraction  of  electrons  that  continue  on  for  all  incident  angles.  These  values  are 
then  multiplied  by  the  electron  flux  at  that  energy  to  provide  as  a  function  of  depth  the  energy 
dose  rate  (Dose  energies),  the  initial  electron  current  (Jo  energies),  and  the  charge  deposition 
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(rho_energies)  at  that  particular  energy.  6  The  loop  is  repeated  for  each  energy  in  the  spectrum. 
After  calculating  for  all  energies,  the  results  are  summed  to  provide  dose  rate  (Dose),  initial 
current  (Jo),  and  deposited  charge  (rho)  as  a  function  of  depth  in  the  dielectric.  These  results  are 
then  provided  to  the  transport  code  shown  in  the  blue  Time  Loop  block  of  Figure  2. 

3.4.2  Input  of  electron  energy  spectrums.  In  order  to  input  the  energy  spectrum,  the  user 
creates  a  comma-delimited  file  (.csv).  The  file  contains  a  minimum  of  two  rows.  The  first  cell 
must  have  always  have  a  zero  entered  in  it  as  a  place  holder  in  the  comma-delimited  file.  After 
that,  each  cell  in  the  first  row  will  list  the  energy  of  a  channel  in  MeV.  The  first  cell  of  the 
second  row  will  give  the  number  of  seconds  that  the  energy  spectrum  being  modeled  is  incident. 
The  rest  of  the  cells  in  the  second  row  will  give  the  electron  flux  in  l/(cm  sr  sec  MeV)  of  the 
energy  channel  listed  above  it.  An  illustrative  example,  with  energy  channels  0.5,  1.0,  2.0,  5.0, 

10,  20  MeV  is  provided  in  Table  1.  In  this  table,  there  are  two  rows  with  differing  electron  flux 
values  that  run  for  100  and  150  seconds.  More  discussion  on  the  input  of  incident  electron  flux 
environments  is  reserved  for  Section  3.5. 

Once  AF-NUMIT2  reads  the  file,  it  calculates  the  size  of  each  energy  channel  in  MeV.  With  the 
exception  of  the  first  and  last  energy  channels,  the  energy  width  is  calculated  from  half  the 
distance  in  MeV  to  the  channel  on  either  side.7  For  the  energy  channels  on  either  end  of  the 
spectrum,  the  half  distance  to  the  next  channel  is  doubled.  In  Table  1  this  width  would  be  0.5 
MeV  for  the  first  channel  and  10  MeV  for  the  last  channel.  The  second  channel  has  a  width  of 
0.75  MeV,  the  third  has  2.0  MeV  and  so  forth.  The  electron  flux  for  a  given  channel  (e.g.,  second 
and  third  lines  in  Table  1)  is  multiplied  by  the  calculated  energy  channel  width  to  provide  the 
electron  flux  of  that  channel  in  l/(cm  sr  sec). 


Table  1:  An  example  of  an  electron  energy  spectrum  input  file 


0 

0.5 

1 

2 

5 

10 

20 

100 

3.67E+04 

4.21E+03 

4.83E+02 

5.55E+01 

6.37E+00 

7.31E-01 

150 

7.34E+04 

8.42E+03 

9.66E+02 

1.11E+02 

1.27E+01 

1.46E+00 

3.4.3  Electron  Beams.  If  the  user  desires  to  model  a  normally  incident  beam  with  a  spectrum  of 
energies,  the  number  of  angular  bins  in  AF-NUMIT2  must  be  set  to  one.  After  that,  the  input  file 
is  prepared  in  almost  the  same  way.  The  difference  is  that  the  electron  flux  for  each  channel  will 
be  in  l/(cm  sec  MeV)  because  the  steradian  units  are  meaningless  for  a  beam.  For  mono- 


6  We  use  parenthesis  to  identify  the  parameter  name  that  appears  in  the  actual  AF-NUMIT2  code. 

7  A  provision  does  exist  in  the  code  to  make  this  calculation  logarithmic.  That  is,  the  width  is  calculated  from  the 
natural  log  of  the  energy  distance  to  the  channels  on  either  side. 
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energetic  beams,  only  one  energy  channel  and  flux  should  be  given  in  the  comma-delimited 
input  file.  In  that  case,  the  electron  beam  flux  will  simply  be  l/(cm  sec). 


3.5  Provisions  for  Time-sequenced  Incident  Electron  Flux  Environments 

AF-NUMIT2  can  easily  handle  a  changing  energy  spectrum  electron  flux  environment.  In  the 
comma-delimited  input  file  additional  rows  are  added  for  new  electron  fluxes.  The  first  cell  in 
each  row  designates  the  time  in  seconds  that  that  flux  will  be  incident.  In  the  example  of  Table  1 
these  values  are  100  seconds  for  the  first  spectrum  and  150  seconds  for  the  second.  The  other 
cells  in  each  row  indicate  the  electron  flux  for  the  corresponding  energy  channel.  In  the  Table  1 
example,  the  second  electron  flux  for  each  energy  channel  happens  to  be  twice  that  of  the  first. 
The  energy  spectrum  loop  will  be  run  for  each  electron  flux  energy  spectrum  (i.e.,  each  row  past 
the  first  row  in  the  energy  spectrum  input  file).  Once  the  dose  rate  and  initial  current  profiles  are 
detennined,  the  transport  code  in  the  time  loop  will  track  the  movement  of  charges  for  the 
designated  length  of  time — in  the  example  100  and  then  150  seconds.  During  this  time,  the 
model  continues  to  simulate  the  input  of  additional  electrons  and  energy  as  specified  in  the  input 
file.  At  the  end  of  the  simulation  time  for  an  electron  flux  energy  spectrum,  the  final  resulting 
electric  field,  charge  distribution,  and  current  profile — as  determined  during  the  time  loop — are 
stored  for  plotting.  The  dose  rate  used  as  determined  by  BeeckenTabFred  is  also  stored.  In  the 
example,  plots  would  be  made  representing  the  situation  at  100  seconds  and  250  seconds. 

If  the  user  does  not  desire  to  change  the  incident  electron  flux  environment,  but  wishes  to  track 
the  development  of  the  calculated  parameters  in  time,  then  in  the  .csv  input  file  the  user  puts  the 
same  flux  in  multiple  rows.  An  updated  plot  will  be  made  of  each  parameter  after  the  time 
chosen  for  each  row  is  completed. 

It  is  important  to  note  what  the  AF-NUMIT2  model  does  when  changing  from  one  electron  flux 
energy  spectrum  to  the  next  (one  row  to  the  next  in  the  input  file).  The  only  carry-over  from  the 
previous  electron  flux  is  the  position  dependent  charge  deposition  profile  (rho).  The  dose  rate 
profile  (Dose)  starts  over  for  the  new  flux,  in  keeping  with  the  thinking  that  radiation  induced 
conductivity  “RIC”  has  no  persistence  or  time  dependence  (see  Section  3.8).  Obviously  a  new 
calculation  must  be  made  of  the  initial  current  (Jo)  because  it  is  due  solely  to  the  incident 
electrons.  The  charge  induced  on  the  electrodes  is  recalculated  on  every  trip  through  the  Time 
Loop  and  is  based  solely  on  the  charge  that  is  deposited  within  the  dielectric.  The  electric  field  is 
also  calculated  each  trip  through  the  Time  Loop  and  it  is  based  solely  on  the  deposited  charge 
and  the  charges  induced  on  the  electrodes. 
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3.6  Extension  of  Model  for  Incident  Electrons  to  Lower  Energy  Range 


A  significant  limitation  of  NUMIT  to  be  addressed  in  the  present  work  was  its  restriction  to 
modeling  only  incident  electrons  with  energies  greater  than  100  keV.  In  the  space  environment 
there  are  normally  significant  fluxes  of  electrons  at  lower  energies,  and  these  are  believed  to 
make  important  contributions  to  the  deep  dielectric  charging  problem.  There  is  nothing  inherent 
in  the  transport  model  part  of  NUMIT  that  excludes  the  use  of  lower  energies.  However,  as 
published,  both  the  energy  deposition  algorithm  EDPOS  [3]  and  the  charge  deposition  model 
FredBell  [4]  are  only  valid  for  incident  electron  energies  0. 1-20  MeV.  The  task  then  is  to 
extend  EDEPOS  and  FredBell  to  lower  energies  or  to  replace  them. 

3.6.1  A  Lower  Incident  Energy  Model.  An  important  effort  has  been  made  by  a  group  at  JPL  to 
address  the  limitations  of  EDEPOS.  [18]  Kim  et  al.  did  Monte  Carlo  electron  transport 
calculations  using  TIGER/ITS3.  They  observed  that  for  a  given  material,  with  incident  energies 
10-100  keV,  the  energy  dose  profiles  all  have  the  same  shape  and  “can  be  reduced  to  a  single 
curve  by  normalizing  them  using  energy  dependent  scaling  factors  for  the  depth  and  the  energy 
deposition.”9 


Naturally,  the  total  of  all  deposited  energy  must  be  dependent  on  the  incident  energy. 
Accordingly,  Tabata  and  Ito  [3]  showed  that  the  integral  of  the  dose  D  is  proportional  to  the 
incident  energy  T0, 


dx  =  fT0 


where/is  a  normalization  factor  that  accounts  for  the  fraction  of  incident  energy  backscattered 
from  the  surface  and  the  energy  deposited  deeper  due  to  the  bremsstrahlung  tail.  Implicit  in  this 
formulation  is  the  use  of  x  as  depth  from  the  surface  of  the  material  written  in  terms  of  g/ctrT. 

Commonly,  the  depth  is  converted  to  normalized  depth  w  =  x/Rex  ,  where  Rex  is  the 
extrapolated  range  of  the  incident  electrons  for  a  given  material  and  energy.  Kim  develops  a 
normalized  dose  profile  g  (w)  and  relates  it  to  the  actual  dose  profile  through  a  vertical  scaling 
factor  Sy\ 


D(x)  =  Sy  g(w ). 


8  Frederickson  does  provide  tables  that  allow  for  interpolation  of  FredBell  lip  to  100  MeV. 

9  Frederickson  suggested  to  this  author  in  2000  that  we  ought  to  be  able  to  make  a  reasonable  approximation  of  the 
dose  profiles  at  lower  energies  by  scaling  down  from  a  higher  energy  using  the  incident  energy  ratio  and  the  ratio  of 
extrapolated  ranges.  Kim’s  detailed  work  validates  this  approximate  approach. 
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Kim  found  an  equation  for  the  nonnalized  dose  profde  g(w ),  but  does  not  provide  the  necessary 
fitting  parameters.  Nevertheless,  we  can  make  progress  because  Kim  has  validated  the  scaling 
principle. 

For  the  specific  case  of  100  keV  incident  energy  electrons  the  EDEPOS  algorithm  claims 
validity.  Thus,  we  choose  D(x)  ->  D100(x ).  In  addition,  because  Kim’s  g(w)  is  valid  for  the 
incident  energy  range  10-100  keV,  we  can  also  write  the  following  expression  for  Diow(x )  at  all 
incident  energies  in  that  range  lower  than  100  keV: 

DiowW  =  Sylow  g(w ). 

By  combining  these  last  two  equations  we  can  eliminate  Kim’s  nonnalized  dose  profile: 

D, „„(*)  =  £>! oo(X>. 


Kim  shows  that 


c  = 
y  R 


fT0 


Therefore,  for  incident  energies  lower  than  Ta=  100  keV,  we  can  also  write: 


yiow 


flow  To low 
R exlow 


By  substitution  we  now  have  an  expression  that  converts  a  dose  profile  at  100  keV  to  a  dose 
profile  at  lower  energies: 


(4) 


Using  their  Monte  Carlo  simulations  over  the  incident  energy  range  10-100  keV,  Kim  finds 
values  for  f  an  equation  for  Rex,  and  a  functional  fit  for  the  normalized  dose  profile  g  (w) . 
Clearly,  this  is  a  significant  step  forward,  and  Kim  demonstrates  how  very  precisely  the 
formulation  fits  some  of  the  Monte  Carlo  simulations.  Unfortunately,  two  fitting  parameters  are 
required  for  Rex  and  six  more  fitting  parameters  are  required  for  g(w).  Although  graphs  of  these 
parameters  were  published,  no  equations  for  the  parameters  were  found.  In  addition,  the 
normalization  factor/ was  presented  as  a  three-dimensional  plot  with  the  proposal  that  a  simple 
fitting  function,  requiring  four  more  fitting  parameters,  may  be  developed  that  could  be  used 
instead  of  their  plotted  f 
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3.6.2  Implementation  in  AF-NUMIT2.  For  the  purposes  of  extending  AF-NUMIT2  to  lower 
energies,  the  work  by  Kim  is  enormously  helpful.  However,  with  a  total  of  twelve  fitting 
parameters  required,  none  of  which  has  an  equation  provided,  it  is  clear  that  the  degree  of 
precision  of  the  work  done  by  Kim  could  not  be  readily  adapted  to  AF-NUMIT2.  Nevertheless, 
the  basic  idea  represented  by  Equation  (4),  derived  from  our  analysis  of  Kim’s  work,  provides  a 
path  to  an  approximation  without  the  necessary  fitting  parameters. 

Tabata  [3]  [19]  claims  a  region  of  validity  for  his  algorithm  for  extrapolated  range  Rex  down  to 
about  0.3  keV.  His  algorithm  for  the  nonnalization  factor  /is  said  to  be  valid  down  to  50  keV. 
Furthennore,  Kim  [18]  reports  that  comparison  of  their/ with  Tabata’s/ down  to  10  keV  yields  a 
deviation  within  ±4  %.  Thus,  using  Tabata’s  algorithms,  we  can  calculate  two  of  the  three  factors 
in  Equation  (4).  The  third  factor  is  a  simple  ratio  of  incident  energies.  All  that  remains  for  an 
approximation  of  the  dose  profile  Dtow(x)  below  100  keV  is  a  good  profile  for  D(x)  at  100  keV. 
Tabata’s  EDEPOS  is  claimed  to  be  valid  down  to  100  keV,  so  that  algorithm  fills  the 
requirement. 

The  BeeckenTabFred  algorithm  for  charge  and  energy  deposition  profiles  branches  at  To  =  100 
keV.  If  the  incident  energy  is  at  or  above  100  keV,  then  a  streamlined,  combined  version  of 
EDEPOS  and  FredBell  is  used.  If,  however,  the  incident  energy  is  below  100  keV,  then  the  low 
energy  branch  is  utilized.  In  this  branch,  the  profiles  at  100  keV  are  calculated  first.  Then, 
Equation  (4)  is  used  to  scale  the  vertical  axis  of  the  Dose,  corresponding  to  the  energy  deposited. 
The  depth  or  horizontal  axis  is  scaled  by  Rexlow/Rexi  which  is  consistent  with  the  normalized 
depth  w  utilized  by  Kim  and  many  others. 

What  about  the  electron  deposition?  The  FredBell  algorithm  provides  a  profile  of  the  incident 
current  as  it  drops  to  zero.  To  find  electron  deposition  from  it,  the  user  simply  takes  the 
derivative.  Therefore,  when  scaling  to  lower  than  100  keV,  there  is  no  reason  to  adjust  the 
vertical  scale.  We  simply  take  the  derivative  and  scale  the  horizontal  axis  in  the  same  manner  as 
for  the  Dose. 
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3.6.3  Results  of  new  Lower  Incident  Energy  Model.  Figure  5  shows  the  results  of  Kim’s  dose 
profiles  determined  by  Monte  Carlo  calculations  (the  only  profiles  provided  using  incident 
energies  below  100  keV).  [18]  In  Figure  6,  the  profiles  were  found  using  the  approximation  of 


Figure  5.  Incident  energy  deposition  plot  made  by  new  JPL  algorithm  and  Monte  Carlo  work 


Dose  Profile:  Carbon 


Figure  6.  Incident  energy  deposition  plots  made  by  BeeckenTabFred  approximate  algorithm 

the  BeeckenTabFred  algorithm.  Although  similar,  the  two  profiles  for  normal  incidence  are  not 
the  same.  The  peaks  are  surprisingly  close  in  magnitude,  but  the  approximate  algorithm  method 
has  the  peaks  forming  slightly  closer  to  the  surface  than  the  Monte  Carlo  calculations.  Why  the 
difference?  There  are  four  possibilities: 
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1 .  The  factor fow  used  in  the  approximate  algorithm  is  Tabata’s,  and  it  is  not  claimed  to  be 
valid  below  T0  =  50  keV. 

2.  The  algorithm  uses  Rex  as  determined  by  Tabata,  which  although  according  to  Tabata  is 
valid  in  the  range  plotted,  is  not  the  same  extrapolated  range  as  used  by  Kim. 

3.  Kim  does  not  provide  the  density  for  the  carbon  used,  and  since  carbon  has  a  range  of 
densities,  the  value  used  in  our  algorithm  was  simply  a  guess. 

4.  The  algorithm  requires  a  good  fit  by  EDEPOS  at  100  keV  because  that  is  the  dose  profile 
that  is  scaled  to  all  the  lower  energies. 

The  first  possibility  would  seem  to  be  an  obvious  problem.  However,  Equation  (4)  only  scales 
the  height  of  the  dose  rate  curve,  not  the  shape.  Therefore,  since  fow  is  only  used  by  the 
BeeckenTabFred  algorithm  in  the  context  of  Equation  (4),  it  cannot  be  the  problem. 

In  the  second  possibility,  the  extrapolated  range  Rex  is  utilized  by  BeeckenTabFred  in  two  ways. 
First,  it  is  used  in  Equation  (4),  which  we  have  already  shown  cannot  be  the  issue.  Second,  it  is 
used  to  scale  the  horizontal  axis  by  Rexiow/Rex-  This  could  be  a  problem,  but  Tabata  claims  Rex 
is  valid  to  energies  way  below  those  considered  here  (~0.3  keV).  Therefore,  it  would  seem 
unlikely  the  discrepancy  comes  from  this  source. 

For  the  third  possibility,  if  different  values  of  density  were  chosen  for  carbon,  the  problem  is 
specific  to  this  case.  The  BeeckenTabFred  algorithm  was  run  for  carbon  with  densities  of  1.8 
g/cm  and  2.26  g/cnr  (and  even  the  unrealistically  extreme  case  of  0.8  g/cm  ).  The  results  were 
identical — not  unsurprising  because  the  Depth  (mg/cm")  indicated  by  the  horizontal  axes  in 
Figure  6  is  scaled  by  the  density. 

The  fourth  and  last  possibility  remains.  Without  any  real  doubt,  the  discrepancy  between 
Tabata’s  EDEPOS  algorithm  at  To  =  100  keV  and  Kim’s  Monte  Carlo  result  is  the  cause  of  the 
differences  observed  between  Figure  5  and  Figure  6.  The  100  keV  line  in  the  lower  plot  comes 
directly  from  EDEPOS  because  that  falls  within  the  validity  range.  Note  it  peaks  at  about  5.5 
mg/cm"  and  the  equivalent  line  in  the  upper  JPL  plot  peaks  at  about  7  g/cm".  When  scaling  the 
horizontal  axis  to  lower  energies,  this  difference  will  be  propagated.  As  discussed  in  the  next 
section,  there  is  reason  to  believe  that  the  JPL  algorithm  by  Kirn  is  better  than  EDEPOS  at  To  = 
100  keV. 


3.7  Verification  of  the  Energy  Deposition  Model  in  NUMIT 

In  order  for  AF-NUM1T2  to  be  useful,  the  depth  profiles  for  the  deposition  of  the  incoming 
electrons  and  their  energies  must  approximate  reality.  As  discussed  earlier,  Tabata  and 
Frederickson  used  limited  experimental  data  and  Monte  Carlo  simulations  to  develop  algorithms 
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that  describe  these  depth  profiles  for  various  materials.  Because  these  algorithms  were  only 
designed  to  be  valid  for  incident  electron  energies  above  100  keV,  the  last  section  was  devoted  to 
describing  the  present  attempt  (the  algorithm  BeeckenTabFred)  to  extend  the  profiles  to  lower 
incident  energies.  But  how  valid  is  this  attempt,  and  for  that  matter,  how  valid  are  the  original 
algorithms? 

In  fact,  why  bother  with  the  approximate  algorithms?  Why  not  run  a  Monte  Carlo  simulation  for 
any  situation  to  be  modeled  by  AF-NUMIT2?  Monte  Carlo  can  be  used  to  track  the  deposition  of 
energy  and  charge  in  the  dielectric  and  from  this  we  can  produce  energy  dose  and  charge 
deposition  depth  profiles.  Unfortunately,  even  on  modem  computers,  Monte  Carlo  simulations 
often  take  hours  to  run  and  require  being  redone  for  every  incident  electron  energy  and/or  change 
of  material.  Such  an  approach  is  impractical  for  modeling  the  simultaneous  incidence  of  a  wide 
energy  spectrum  of  electrons — the  reason  NUMIT  is  being  redesigned.  More  importantly,  AF- 
NUMIT2  is  a  tool  for  investigating  various  factors — materials,  conductivities,  temperatures,  flux 
spectra,  and  combinations  of  flux  spectra — in  an  effort  to  determine  what  conditions  make 
electrostatic  discharge  likely.  There  are  so  many  possible  combinations  of  parameters  that  slow 
processing  would  make  such  investigations  unwieldy. 

Because  the  BeeckenTabFred  algorithm  is  based  on  prior  work,  we  believe  it  should  be 
sufficiently  accurate,  but  the  only  known  basis  of  comparison  is  one  graph  published  by  JPL  [18] 
for  energy  deposition  in  the  material  carbon  (reproduced  in  Figure  5).  The  goal  here  then  is  to 
produce  Monte  Carlo  simulations  that  will  be  able  to  determine  the  efficacy  of  the  extended 
algorithm.  Hopefully,  this  simple,  fast  algorithm  will  at  least  approximate  the  energy  and  charge 
deposition  profiles  in  various  materials  at  lower  incident  electron  energies. 

3.7.1  Setup  of  the  Monte  Carlo  Method.  MCNP6  Beta  Version  2  was  used  to  run  Monte  Carlo 
(MC)  simulations  to  find  the  energy  and  charge  deposition  profiles  in  different  materials  with 
incident  electron  energies  between  10  and  100  keV.  In  order  to  run  a  simulation  using  MCNP6, 
input  files  need  to  be  created  with  the  correct  problem  and  geometry  setup.  The  material  is 
divided  into  many  surfaces,  all  stacked  up  in  series  with  a  small  “segment”  in  between  each 
surface.  MCNP6  fires  electrons  at  the  resulting  mesh.  The  segment  width  can  be  varied  as 
necessary  to  produce  the  desired  resolution  of  the  energy  curves.  More  surfaces  with  smaller 
segments  in  between  construct  a  finer  mesh  which  will  in  turn  result  in  greater  resolution — 
important  particularly  for  the  sub  50  keV  energies.  However,  more  surfaces  also  greatly 
increases  the  computation  time  for  the  simulation  so  that  a  coarser  mesh  must  be  used  at  higher 
energies  where  greater  penetration  depths  are  modeled.  The  cell  is  comprised  of  a  single 
material  which  is  defined  by  the  effective  atomic  number,  the  effective  mass  number,  and  the 
material  density.  The  total  energy  per  gram  deposited  in  each  segment  is  tallied  by  MCNP6  and 
written  to  an  output  file. 
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After  the  geometry  of  the  problem  is  set  up,  the  electron  source  needs  to  be  defined.  For  these 
MC  simulations  we  sent  one  million  electrons  into  the  material  and  MCNP6  tracked  their 
movement  and  the  deposition  of  their  energy  within  the  mesh.  Once  an  electron  deposits  its 
energy  and  is  left  with  100  eV  or  less,  the  tracking  stops.  Monte  Carlo  simulations  were  run  for  a 
few  different  materials,  calculating  both  the  energy  and  charge  deposition.  They  are  described 
below. 


3.7.2  Examination  of  Low  Energy  Algorithm  for  Energy  Deposition.  Monte  Carlo 
simulations  of  energy  deposition  have  been  published  for  carbon  by  JPL  and  were  reproduced  in 
Figure  5  along  with  our  extended  BeeckenTabFred  algorithm  results  in  Figure  6.  As  a  baseline, 
we  ran  MNCP  on  the  same  situation.  Figure  7  is  the  result  of  MC  simulations  produced  by  David 
Dixon,  a  graduate  student  at  the  University  of  New  Mexico  who  introduced  us  to  MCNP.  Using 
what  we  learned  from  Dixon,  we  were  able  to  produce  a  similar  result,  Figure  8  for  carbon  (Z=6, 
A=12,  p=2.1  g/cnr).  Comparison  of  Figures  5,  7,  and  8  shows  the  similarity  between  the  three 
Monte  Carlo  plots.  The  peaks  are  all  in  the  same  place,  although  of  slightly  different  heights.  The 
differing  heights  are  apparently  due  to  subtle  differences  between  versions  of  MNCP6  and  the 
TIGER/ITS3  used  by  Kim  at  JPL.  We  used  beta  version  2  on  our  plot,  but  we  have  also  used 
identical  inputs  to  compare  the  results  to  beta  version  3  and  found  the  resulting  vertical  scales  to 
be  different.  We  are  not  sure  what  version  David  Dixon  used  for  his  plot,  but  we  believe  it  was 
MCNP5.  Regardless,  the  main  purpose  here  is  to  verify  the  shape  of  the  dose  curve. 


As  mentioned  at  the  end  of  Section  3.6,  the  peaks  of  the  dose  in  the  extended  algorithm  plots  of 
Figure  6  are  about  20  %  closer  to  the  surface  than  the  MC  plots.  Since  the  MC  plots  of  Figures  5, 


Dose  Profile:  Carbon  (Dixon) 
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Figure  7.  Monte  Carlo  simulations  of  normal  incidence  on  Carbon  by  David  Dixon 
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7,  and  8  are  all  consistent,  it  looks  bad  for  our  extended  algorithm.  However,  the  difference  holds 
true  for  the  100  keV  case  which  is  not  part  of  the  extended  algorithm  but  a  result  of  the  original 
EDEPOS  algorithm.  Thus,  it  seems  that  the  EDEPOS  algorithm  at  100  keV  is  not  as  accurate  as 
Kim’s  algorithm  which  he  fit  to  their  MC  simulations  at  100  keV.  It  would  be  interesting  to 


Dose  Profile:  Carbon  (MC) 


Figure  8.  Monte  Carlo  simulations  of  normal  incidence  on  Carbon 


follow  this  difference  between  modern  MC  results  and  EDEPOS  through  to  higher  energies  and 
in  different  materials.  Time,  however,  is  a  limiting  factor  and  EDEPOS  has  been  thought  for 
decades  to  be  sufficiently  accurate  for  modeling  purposes. 


Teflon  was  the  second  material  that  was  used  to  compare  MC  simulations  with  the  new 
BeeckenTabFred  algorithm.  Energy  deposition  curves  were  produced  with  MCNP  for  incident 
electron  energies  from  10  to  100  keV  in  increments  of  10  keV,  see  Figure  9.  In  this  case,  the 
same  segment  width  of  1 .0  pm  was  used  in  the  geometry  setup  for  most  energies,  although 
smaller  segments  were  needed  for  the  lower  energies  of  10-40  keV. 
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Dose  Profile: Teflon  (MC) 


Figure  9.  Monte  Carlo  simulations  of  normal  incidence  on  Teflon 


The  BeeckenTabFred  algorithm  was  then  run  for  Teflon  using  the  effective  atomic  number  8.28 
and  effective  atomic  weight  17.25  with  mass  density  2.2  g/cm3,  see  Figure  10.  The  results  shown 
in  Figure  9  and  Figure  10  are  quite  similar,  but  in  order  to  get  an  accurate  picture  of  the 
differences,  individual  energy  curves  must  be  plotted  on  the  same  graph  as  shown  in  Figure  1 1 

The  BeeckenTabFred  algorithm,  which  for  this  energy  identical  to  EDEPOS,  is  slightly  offset  to 
towards  the  surface  compared  to  the  MC  results. 10  It  is  clear  that  the  EDEPOS  algorithm  predicts 
a  slightly  shallower  deposition  of  energy  at  100  keV  than  our  current  MC  simulations  show. 
Because  BeeckenTabFred  is  a  scaling  of  the  100  keV  case,  this  slight  offset  propagates  to  the 
rest  of  the  energies  that  were  compared  for  Teflon  and  carbon. 

Figure  12  provides  further  comparisons  between  MC  simulations  and  the  BeeckenTabFred 
algorithm  that  extends  EDEPOS  to  lower  energies.  Although  the  mesh  size  for  the  20  keV  case 
was  apparently  too  coarse  for  good  results,  it  appears  again  that  the  MC  simulations  are 
indicating  slightly  deeper  energy  deposition  compared  to  the  algorithm — consistent  with 
expectations  based  on  the  difference  at  100  keV. 


10  The  results  from  the  MC  simulation  were  vertically  scaled  by  a  factor  of  1.47  in  order  to  more  effectively  compare 
the  two  deposition  curves.  Vertical  scales  are  difficult  to  correlate  in  this  type  of  modeling.  At  this  point  in  the 
research,  the  shape  of  the  distribution  is  of  primary  concern.  Furthermore,  we  have  observed  that  performing  exactly 
the  same  simulation  with  beta  versions  2  and  3  of  MCNP6  resulted  in  different  vertical  scales. 
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Energy  Deposition:  Teflon  (Beecken) 


Figure  10.  Results  of  BeeckenTabFred  algorithm  for  normal  incidence  on  Teflon 


1 00  keV 


Figure  11.  Comparison  of  BeeckenTabFred  algorithm  with  MC  simulation 

The  slight  ripple  of  the  MC  curve  in  Figure  1 1  for  100  keV  was  investigated  in  order  to 
determine  whether  or  not  the  mesh  size  was  the  issue.  The  “stair  step”  nature  of  the  graph  would 
seem  likely  to  be  caused  by  simulating  a  continuous  process  with  a  discrete  algorithm,  therefore 
the  MC  simulation  was  repeated  for  Teflon  with  100  keV  incident  electron  energies  but  with  the 
mesh  size  cut  in  half  from  the  previous  simulation  run.  The  two  curves  have  the  same  ripple,  and 
the  smaller  mesh  size  did  not  improve  the  resolution  of  the  curve.  Perhaps  an  even  smaller  mesh 
would  succeed  in  eliminating  this  effect,  but  this  seems  unlikely.  Another  possible  solution 
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Figure  12.  Detailed  comparisons  of  algorithm  with  MC  at  50  and  20  keV. 


would  be  to  run  the  MC  simulation  again  with  many  more  incident  electrons.  Similar  ripples 
show  up  in  Dixon’s  independent  results  in  Figure  7. 

Comparisons  between  Monte  Carlo  simulations  and  the  BeeckenTabFred  algorithm  were  also 
made  for  aluminum  which  has  greater  atomic  number  (13),  atomic  weight  (27),  and  density  (2.7 
g/cm  ).  The  results  showed  the  same  trend  identified  above. 

3.7.3  Examination  of  Low  Energy  Algorithm  for  Charge  Deposition.  Although  it  was  not  part 
of  the  statement  of  work,  it  is  important  to  also  investigate  the  validity  of  the  BeeckenTabFred 
algorithm  in  terms  of  charge  deposition.  Monte  Carlo  simulations  were  performed  at  energies 
between  10  keV  and  100  keV  for  carbon,  Teflon,  and  aluminum.  For  aluminum,  we  looked  only 
at  where  the  electrons  stopped  due  to  the  material  stopping  effect  and  inferred  a  current  from 
that.  Effectively,  this  makes  the  aluminum  non-conducting  for  the  purposes  of  the  simulation,  but 
aluminum  is  such  a  standard  reference  material  for  radiation  shielding  it  was  an  appropriate 
choice. 

Figure  13  illustrates  the  results  for  aluminum  (non-conducting)  at  100,  50,  and  10  keV.  As 
expected  from  the  energy  deposition  results,  the  electron  depositions  from  the  algorithm  were 
slightly  closer  to  the  surface  for  100  and  50  keV.  Surprisingly,  the  10  keV  case  had  MC 
depositing  the  electrons  closer  to  the  surface.  The  result  is  interesting  because  for  carbon  and 
Teflon  the  trend  had  held.  Probably,  what  is  really  being  indicated  by  these  differences  is  that  the 
approximations  made  by  the  algorithm  worsen  as  the  extreme  low  energy  of  10  keV  is 
approached. 


ii 


The  results  of  these  comparisons  are  documented  more  thoroughly  in  the  21  Month  After  Award  report. 
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Remaining  Fraction  of  Incoming  Electrons  Remaining  Fraction  of  Incoming  Electrons  Remaining  Fraction  of  Incoming  Electrons 


Primary  Electron  Current:  Aluminum  100  keV 
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Primary  Electron  Current:  Aluminum  50  keV 


Depth  (mm) 


Primary  Electron  Current:  Aluminum  1 0  keV 


Figure  13.  Comparisons  of  electron  deposition  results:  BeeckenTabFred  algorithm  vs.  Monte  Carlo 
simulations 
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3.7.4  Examination  of  Algorithms  for  Isotropic  Incidence.  The  goal  in  redesigning  NUMIT 
into  a  tool  useful  for  modeling  the  spacecraft  environment  necessarily  included  accounting  for 
electron  incidence  at  various  angles  other  than  perpendicular  to  the  surface.  The  method  chosen 
was  described  in  Subsection  3.3.2.  In  an  effort  to  test  the  validity  of  the  method,  we  ran  MC 
simulations  for  a  large  number  of  circumstances.  This  task  is  rather  involved  because  there  are  so 
many  factors.  We  used  MCNP6  Beta  Version  3  for  incident  energies  of  10,  20,  50,  100,  200, 

500,  and  1000  keV  on  carbon  and  aluminum.  Angles  used  were  0°,  30°,  and  60°.  Results  were 
found  for  both  energy  and  charge  depositions.  All  totaled,  for  all  of  these  possible  combinations, 
there  would  be  84  plots.  We  ran  82  of  them,  requiring  hundreds  of  hours  of  computation  time. 

The  sheer  volume  of  results  was  somewhat  overwhelming,  but  some  general  observations  can  be 
made.  Many  of  the  MC  simulations  were  very  similar  to  the  algorithm.  Normally,  the  peaks  were 
relatively  close.  The  similarity  deteriorated  as  the  angle  from  the  normal  was  increased. 

However,  due  to  enhanced  backscatter  at  greater  angles  there  are  fewer  electrons  entering  the 
material  at  these  angles.  Therefore  the  discrepancy  is  less  of  a  problem. 

Unfortunately,  this  effort  at  validation  was  very  time  consuming  and  not  really  conclusive.  It  did 
seem  to  indicate  strongly  that  the  algorithm  approach  being  taken  is  valid  as  a  first  order 
approximation.  In  retrospect,  however,  some  things  should  have  been  done  differently.  First, 
given  the  difference  in  energy  deposition  quantities  between  MCNP6  Beta  Version  2  and 
Version  3,  it  would  have  been  better  to  use  MCNP5  or  wait  until  MCNP6  was  finalized. 
Secondly,  the  algorithm  runs  that  were  plotted  as  comparisons  on  the  same  graphs  were  not  done 
with  backscatter  accounted  for  in  the  algorithm.  Since  MCNP  does  take  backscatter  into 
consideration,  the  comparisons  on  the  vertical  scales  (energy  or  number  of  electrons  deposited) 
were  compromised  from  the  start.  This  inherent  difference  in  the  basis  of  comparison  will  be 
more  of  an  issue  at  greater  incident  angles  which  have  more  backscatter — precisely  where  the 
discrepancy  increased. 

Further  consideration  of  the  problem  has  resulted  in  the  conclusion  that  the  accuracy  of  the 
algorithm  approximation  for  a  single  angle  is  not  only  inconsistent,  but  it  is  really  not  the 
important  comparison.  Since  it  is  isotropic  fluxes  that  are  being  modeled,  the  only  thing  that 
really  counts  is  how  the  result — after  summing  through  all  the  possible  angles — compares  with 
reality.  Thus,  the  appropriate  approach  is  to  do  Monte  Carlo  simulations  with  isotropic  incidence 
of  electrons  instead  of  using  a  beam  set  at  various  angles. 

Some  preliminary  Monte  Carlo  studies  of  isotropically  incident  electrons  have  been  done  in 
collaboration  with  David  Barton  of  AFRL/RVB.  Barton  used  MCNP5  to  simulate  isotropic 
electrons  incident  on  aluminum  at  various  energies.  These  simulations  included  backscatter  of 
electrons  off  the  surface.  The  electron  deposition  profiles  for  five  different  incident  energies  are 
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plotted  as  points  in  Figure  14.  The  solid  lines  in  the  figure  are  the  results  of  the  BeeckenTabFred 
algorithm  (with  the  backscatter  feature  on)  used  in  AF-NUMIT2. 

Barton  also  found  the  energy  deposition  profiles  in  aluminum  with  MCNP5  for  the  same  incident 
energies  of  isotropic  fluxes.  The  results  of  both  the  algorithm  and  the  MC  simulations  were  put 
on  the  same  scale  with  the  same  units  and  plotted  in  Figure  15.  The  points  are  MC  simulation 
data  and  the  solid  lines  are  the  AF-NUMIT2  deposition  algorithm. 

Although  it  might  appear  that  the  correlation  between  the  MC  simulations  and  the  algorithm  is 
not  impressive,  never  before  have  MC  simulations  been  used  to  find  depth  profiles  of  charge  and 
energy  deposition  for  isotropically  incident  electrons  of  several  energies.  Furthermore,  no 
algorithm  besides  BeeckenTabFred  is  known  to  exist  that  estimates  deposition  depth  profiles 
from  isotropic  incidence,  much  less  for  a  large  range  of  materials  and  energies.  This  algorithm 
was  developed  without  the  aid  of  MC  simulation  or  lab  data — it  is  the  result  of  an  analytical 


Figure  14.  Monte  Carlo  simulations  of  electron  deposition  vs.  BeeckenTabFred  algorithm 
for  isotropically  incident  elctrons  on  aluminum 

extension  of  a  normal  incidence  algorithm.  The  fact  that  these  two  independent  results  have  a 
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Figure  15.  Monte  Carlo  simulations  vs.  BeeckenTabFred  algorithm  for 
isotropically  incident  electrons 


worst  case  discrepancy  that  only  in  some  places  approaches  50  %  is  remarkable.  Given  the  other 
uncertainties  involved  with  deep-dielectric  charging  in  spacecraft,  it  would  seem  that  the 
algorithm  is  an  appropriate  approximation.  At  the  same  time,  these  results  provide  a  clear  path 
forward.  Monte  Carlo  simulations  of  isotropic  incidence  need  to  be  made  for  a  wide  range  of 
energies  on  a  large  selection  of  materials.  Then,  a  better  algorithm  might  possibly  be  developed 
to  insert  into  AF-NUMIT2. 


3.8  Development  of  Time-Delayed  and  Temperature-Dependent  Radiation  Induced 
Conductivity  Algorithms 

Development  of  time-dependent  RIC  and  temperature-dependent  RIC  algorithms  were  originally 
thought  to  be  an  important  part  of  this  work.  However,  during  a  25  July  2012  meeting  at 
AFRL/RVBX,  these  parts  of  the  original  proposal  were  reconsidered.  Some  at  the  meeting 
characterized  such  algorithms  as  relatively  unimportant  compared  to  the  other  goals.  Another 
participant  raised  the  possibility  that  including  such  features  might  be  dangerously  misleading. 
Users  might  come  to  believe  that  AF-NUMIT2  is  correctly  accounting  for  temperature  and  time 
dependence,  when  at  best  the  simulation  could  not  be  more  than  a  guess. 
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In  the  case  of  temperature  dependence,  it  is  very  clear  that  there  can  be  a  significant  effect. 
Ferguson  et  al.  [20]  have  shown  experimentally  that  conductivities  can  change  by  orders  of 
magnitude  with  temperature.  However,  the  lack  of  sufficient  data  over  a  wide  enough  range  for 
important  materials  is  a  major  obstacle  to  creating  an  algorithm.  As  regards  time  dependence, 
although  there  exist  clear  possible  approaches  to  an  algorithm  that  includes  this  effect,  it  is 
almost  impossible  to  find  the  time  constants  required  for  such  models. 

After  much  discussion,  the  consensus  was  that  if  AF-NUMIT2  users  wish  to  incorporate 
temperature  affects,  this  would  be  best  done  by  simply  having  the  users  input  the  RIC 
coefficients  and  dark  coefficients  that  they  believe  exist  at  the  temperatures  they  desire  to  model. 
Changing  those  coefficients  as  the  simulation  progresses  is  rather  easily  accomplished,  and  that 
would  give  the  users  complete  control  of  their  temperature  modeling. 


4.0  RESULTS  AND  DISCUSSION 
4.1  Preliminary  Simulation  Results 

In  order  to  test  AF-NUMIT2,  we  obtained  preliminary  results  by  running  it  extensively  with 
realistic  electron  flux  energy  spectra.  Figure  16  shows  various  spectra  obtained  by  CRRES  that 
were  downloaded  and  formatted  into  input  files  suitable  for  the  simulation.  These  data  were 
used  in  a  variety  of  combinations  as  inputs  to  AF-NUMIT2  in  order  to  predict  the  maximum 
electric  field  within  Kapton  after  at  least  8  days  of  exposure — sufficient  time  in  almost  all  cases 
for  equilibrium  to  be  achieved.  The  simulations  normally  took  about  two  minutes  to  run  for  each 
simulated  day. 


12  The  author  is  grateful  to  W.  Robert  Johnston  of  AFRL/RVBX  for  assistance.  The  data  were  obtained  in  two 
overlapping  sets  during  each  minute  of  collection.  These  data  are  available  at  ftp://virbo.org/users/johnston/crres. 
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Figure  16.  CRRFS  data  used  to  run  simulations  in  AF-NUMIT2 

1  T 

A  wide  range  of  possible  values  exist  for  the  dark  conductivity  and  RIC  coefficients  in  Kapton.  ' 
Dark  conductivities  might  be  as  high  as  6><10’  or  as  low  as  4x10’"  1/(Q  cm).  The  possible  RIC 
coefficients  have  almost  as  large  a  possible  range,  from  5x10’' 6  to  as  low  as  5X 10’ 18  (sec/(D  cm 
rad)).  Figure  17  shows  the  results  of  a  simulation  with  AF-NUMIT2.  In  this  simulation,  four 
days  of  the  “hard”  spectrum  from  Figure  16  were  run,  followed  by  four  days  of  “soft”  spectrum. 
For  all  8  days  we  used  a  dark  conductivity  of  7X 10'  l/(f2  cm)  and  a  RIC  coefficient  of  4.5 x  10’ 

(scc/(Q  cm  rad)).  Each  day  is  plotted  as  a  separate  line  in  these  graphs.  Because  equilibrium  was 
essentially  achieved  for  the  “hard”  spectrum  after  one  day,  all  four  days  are  plotted  on  top  of 
each  other.  Then  the  spectrum  shifted  to  “soft”  for  the  remaining  for  plots  corresponding  to  the 
situation  at  the  end  of  each  day’s  simulation.  The  dose  rate,  or  energy  deposition  curves,  clearly 
reflect  these  two  different  input  spectrums.  The  electron  current  profile  is  always  flat,  indicating 


13  Data  provided  in  private  communication  with  JR  Dennison  of  Utah  State  University. 
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equilibrium  was  achieved  after  one  day  for  both  spectra.  The  maximum  electric  field  did  not 
approach  the  typical  risk  level  for  discharging  (105  V/cm)  before  equilibrium  was  reached  at  the 
end  of  the  first  day,  at  the  end  of  four  days,  or  at  the  end  of  the  total  eight  day  simulation.  In  fact, 
based  on  this  simulation,  there  appears  to  be  no  chance  of  discharge-level  electric  fields  being 
reached. 


Both  of  the  conductivity  values  in  this  simulation  were  at  the  high  end.  We  would  expect  high 
conductivity  to  lead  to  larger  bleed-off  of  deposited  charge.  Therefore,  it  becomes  appropriate  to 
explore  the  entire  range  of  conductivites.  Over  thirty  simulation  runs  were  made  with  various 
combinations  of  incident  energy  spectra,  dark  conductivities,  and  RIC  conductivity  coefficients. 
Displayed  in  Figure  18  are  the  results  of  a  simulation  run  made  with  the  same  hard-soft  energy 
spectra  over  the  same  number  of  days.  This  time,  however,  we  chose  to  decrease  the  RIC 
coefficient  by  two  orders-of-magnitude  to  5x  10'  (scc/(Q  cm  rad)),  at  the  low  end  of  the 
possible  value  range.  The  dark  conductivity  was  also  increased  by  two  orders-of-magnitude  to 
7x  10'“  \/(Q  cm),  but  there  is  reason  to  believe  that  this  did  not  affect  the  results.  During  the 
eight  days,  the  simulation  did  not  reach  equilibrium.  Instead  the  electric  field  increased  until  it 
reached  approximately  4.4x1 05  V/cm,  beyond  the  electric  field  strength  commonly  understood  to 
cause  discharge. 
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Figure  18.  Simulation  results  with  lower  RIC  after  4  days  of  "hard"  and  4  days  of  "soft" 


In  a  report  of  reasonable  length,  space  does  not  permit  the  display  of  all  the  AF-NUMIT2 
simulation  runs  that  have  been  made  with  various  paramters.  Perhaps  more  helpful  anyway  is  the 
summary  chart  in  Figure  19  of  the  maximum  electric  field  achieved  during  these  runs.  All  runs 
were  done  using  the  CRRES  spectra  from  Figure  16.  The  solid  markers  correspond  to  “hard- 
soft”  or  just  “soft”  spectra  runs.  The  hollow  markers  were  simulation  runs  made  using  a  spectra 
for  the  entire  run  that  was  not  designated  as  “hard”  or  “soft.”  Most  runs  were  for  8  days  total, 
with  a  few  for  10  days.  The  runs  were  long  enough  for  equilibrium  to  be  effectively  achieved  or 
until  the  electric  field  strength  exceeded  the  assumed  electrostatic  threshold  strength  of  105 
V/cm. 

The  results  showed  a  dramatic  and  consistent  dependence  on  the  value  of  the  RIC  coefficient, 
increasing  very  proportionately  to  the  reciprocal  of  the  coefficient.  As  stated  earlier,  one  should 
expect  that  less  RIC  would  cause  greater  charging.  Lower  conductivity  will  result  in  less  charge 
bleeding  off  the  dielectric.  The  result  is  a  higher  electric  field.  Two  runs  (indicated  with  the  solid 
blue  square  markers)  were  made  with  unrealistically  high  dark  conductivity  in  order  to  verify  this 
understanding.  The  higher  dark  current  clearly  caused  charge  bleed-off  that  negated  the  effect  of 
the  low  RIC  conductivity — these  are  the  only  points  that  do  not  fall  on  the  straight  trendline. 
Despite  the  expectation  that  higher  fields  will  exist  with  lower  conductivity,  the  consistency  with 
which  this  relationship  is  observed,  regardless  of  flux  spectra,  seems  stunning  to  this  author. 
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Figure  19.  Summary  of  AF-NUMIT2  realistic  simulation  results 

After  these  runs  were  performed,  additional  checks  of  the  simulation  time  step  size  and  spatial 
depth  resolution  were  made.  Using  a  smaller  time  step  made  no  difference  in  the  results 
whatsoever — strong  indication  that  an  appropriate  choice  had  been  made.  However,  particularly 
at  lower  incident  energies,  slightly  different  results  were  obtained  when  the  model  used  a  finer 
spatial  depth  resolution.  Apparently,  greater  spatial  resolution  (moving  from  a  spatial  bin  size 
from  10  pm  to  1  pm)  allowed  the  modeling  of  more  charge  closer  to  the  surface.  As  a  result,  this 
charge  was  able  to  bleed  off  more  rapidly  causing  a  reduction  in  the  maximum  electric  field  by 
as  much  as  20  %.  Unfortunately,  increasing  the  resolution  proportionately  increases  the  actual 
run  time  so  that  the  normal  two  minutes  of  run  per  simulated  day  becomes  twenty  minutes.  More 
work  needs  to  be  done  to  find  the  optimal  tradeoff  between  spatial  resolution  and  practical  run 
times — perhaps  the  time  step  size  could  be  increased  to  compensate  for  greater  resolution. 
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4.2  Other  Issues  Discussed  with  AFRL/RVBX  but  not  in  the  SOW 


4.2.1  Secondary  Electron  Yield.  When  high  energy  electrons  strike  the  surface  of  a  dielectric 
they  can  penetrate  or  backscatter,  but  they  can  also  eject  electrons  that  were  part  of  the  dielectric. 
The  ejection  of  secondary  electrons  from  near  the  surface  of  the  dielectric  certainly  could  cause  a 
positive  charging  of  the  surface.  Such  an  effect  could  be  important  in  explaining  or  predicting  the 
buildup  of  an  electric  field — particularly  near  the  surface  and  particularly  at  very  low  incident 
electron  energies  (<  10  keV).  Unfortunately,  to  the  author’s  knowledge,  not  enough  is  known 
about  secondary  electron  yield  to  be  incorporated  into  a  useful  model.  In  addition,  such  a  surface 
effect  is  not  really  part  of  a  deep-dielectric  charging  model. 

4.2.2  The  Effect  of  Internal  Electric  Fields  on  Incident  Electrons.  Neither  the  original 
NUMIT  nor  AF-NUMIT2  incorporates  in  the  model  any  retarding  effect  a  growing  electric  field 
within  the  dielectric  may  have  on  later  incident  electrons.  The  depth  profiles  of  charge  and 
energy  deposition  ignore  the  existence  of  internal  fields.  Many  discussions  have  occurred 
regarding  the  validity  of  this  approximation. 

In  an  earlier  report,14  this  author  presented  an  order-of-magnitude  calculation  that  indicated  any 
effect  of  an  internal  electric  field  in  NUMIT  would  be  less  than  10  %  until  the  threshold  for 
electrostatic  discharge  is  reached.  However,  in  that  same  discussion  it  was  pointed  out  that 
electrons  at  lower  incident  energies  than  NUMIT  was  modeling  might  be  affected  more  greatly. 
Part  of  the  present  work  is  to  extend  NUMIT  to  lower  energies,  therefore  this  shortcoming  of  the 
NUMIT  approach  may  now  be  more  significant. 

The  significance  of  internal  fields  at  lower  incident  energies,  however,  is  not  universally 
acknowledged.  For  example,  Frederickson,  Woolf,  and  Garth,  wrote:  “For  electrons  well  below 
100  keV,  where  the  material  stopping  power  dominates,  the  field  effects  are  less  important.  It  has 
been  observed  for  electrons  above  1  MeV  that  their  trajectories  are  foreshortened  in  charged 
insulators.’’  [21]  Regardless  of  whether  the  issue  of  internal  fields  is  more  prominent  at  low  or 
high  incident  electron  energies,  incorporating  it  into  AF-NUMIT2  was  not  part  of  the  proposed 
work. 

4.2.3  A  Better  Model  of  Charge  and  Energy  Deposition  Profiles.  As  discussed  in 
Subsection  3.7.4,  it  is  clear  that  the  best  approach  for  detennining  the  charge  and  energy 
deposition  profiles  is  to  do  Monte  Carlo  simulations  for  isotropically  incident  electrons  rather 
than  piecing  together  simulations  at  individual  angles  of  incidence.  This  approach  would  also 
incorporate  backscatter 


14  Brian  Beecken,  “Development  of  the  NUMIT  Simulation  for  Modeling  Deep-dielectric  Charging  in  the  Space 
Environment,”  presented  to  AFRL/RVBX,  August  13,  2009. 
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off  of  the  material’s  surface.  It  does  require,  however,  a  totally  new  algorithm  that  would  fit  to 
many,  many  MC  simulations,  thereby  completely  replacing  EDEPOS  and  FredBell  rather  than 
just  modifying  them  as  done  here. 

David  Barton  of  RVBX  has  begun  a  series  of  isotropic  MC  simulations  at  a  wide  range  of 
energies  for  a  wide  range  of  materials.  The  hope  is  to  eventually  compile  enough  MC  results  that 
an  entirely  new  algorithm  can  be  devised  and  then  implemented.  Creating  an  algorithm  that  fits 
MC  energy  and  charge  deposition  profiles  for  a  variety  of  materials  at  a  variety  of  energies 
would  bypass  the  need  to  assume  things  like  the  isotropy  of  dielectrics  and  the  need  to  sum 
multiple  angular  bins.  In  addition,  since  backscatter  is  included  in  the  MC  simulations,  it  would 
automatically  be  included  in  the  algorithm  that  is  fit  to  those  simulations.  Such  an  approach 
would  create  the  basis  of  AF-NUMIT3  and  be  an  exciting  step  forward  in  modeling  of  the  deep- 
dielectric  charging  of  spacecraft. 


5.0  CONCLUSIONS 

The  work  on  development  of  AF-NUMIT2  has  been  performed  and  completed  as  proposed.  The 
only  exception  is  the  inclusion  in  the  model  of  time-dependent  and  temperature-dependent 
radiation  induced  conductivity.  These  dependencies  were  deemed  by  AFRL/RVBX  as  of  lesser 
importance  and  their  incorporation  could  be  a  negative  feature  due  to  the  potential  of  misleading 
future  users. 

Going  beyond  the  proposed  work,  AF-NUMIT2  has  been  tested  extensively  on  Kapton  over  a 
wide  range  of  possible  conductivities,  both  dark  and  RIC,  with  the  use  of  numerous  electron 
energy  flux  spectra  that  are  taken  directly  from  the  data  compiled  from  CRRES  measurements. 
An  inverse  proportional  relationship  between  the  RIC  conductivity  and  the  maximum  internal 
electric  field  has  been  clearly,  and  rather  unexpectedly,  observed  during  the  testing. 
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APPENDIX  A  -  User’s  Guide 


This  is  a  quick  start-up  guide  to  using  AF-NUMIT2.  All  inputs  that  it  refers  to  can  be  easily 
found  in  the  first  100  lines  of  the  code.  For  convenience,  that  code  is  reproduced  in  Appendix  B. 

Getting  Started 

AFNUMIT2.m  is  the  main  AF-NUMIT2  program.  Open  and  run  this  in  Matlab,  but  first  verify 
the  following. 

4  Program  Function  Files  that  must  be  in  the  same  folder  as  AFNUMIT2.m: 

1.  BeeckenTabFred.m  Function  predicts  electron  &  energy  deposition  as  a  function  of  angle 

2.  getparameterB.m  Function  finds  value  of  parameter  B  for  electron  deposition 

3.  backscatter.m  Function  that  guesses  at  backscatter  of  incident  electrons 

4.  num  der  .m  Function  that  takes  a  numerical  derivative 

Once  the  user  verifies  that  the  main  program  and  its  four  function  files  are  all  in  the  same  folder, 
then  the  AFNUMIT2.m  file  must  be  opened.  The  first  choice  that  must  be  made  is  whether  the 
user  is  starting  a  new  run,  or  wants  to  add  additional  simulation  time  to  a  run  that  was  just 
completed.  The  default  is  a  new  run,  and  so  the  first  line  of  code  is  [more  =  'n'].  Here  the 
choice  is  set  to  “n”  which  stands  for  new  or  no.  If  the  simulation  run  is  to  pick  up  where  a 
previous  run  left  off,  the  “n”  is  replaced  with  a  “y”  for  yes. 

Inputs  Required  From  User 

The  parameters  listed  below  must  be  typed  into  the  code.  If  they  are  not  typed  into  the  code,  the 
last  used  parameters  will  automatically  be  reused.  The  obvious  advantage  of  this  approach  is  that 
the  user  can  easily  make  succeeding  runs  by  simply  using  the  file  again,  only  changing  what  is 
needed.  Plus,  it  will  always  be  clear  what  parameters  where  chosen  on  the  run  just  completed. 

Required  Characteristics  of  Dielectric  Material 

Thickness  (cm) 

Effective  atomic  number 
Effective  atomic  weight 
Density  (g/cm  ) 

Dark  Conductivity  [l/(f2  cm)] 

Coefficient  of  RIC  [sec/(Q  cm  rad)] 

Permittivity  (F/m) 
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Resolution  Control  Settings 

1 .  Isotropic  vs.  Beam:  Either  the  incident  electrons  will  be  modeled  as  arriving  isotropically,  or 
they  will  arrive  in  a  beam.  If  the  choice  is  isotropic,  then  the  number  of  angular  bins  must  be 
set,  thereby  detennining  the  angular  resolution  for  isotropic  incidence.  Obviously  choosing 

a  greater  number  will  also  slow  the  simulation  run  time.  It  seems  that  90  bins  works  well.  If 
however,  the  choice  is  a  beam  to  simulate  a  lab  setting,  then  the  number  of  angular  bins 

must  be  set  to  1.  If  it  is  a  mono-energetic  beam,  then  the  flux  units  must  be  l/(cm~  sec),  but 

2 

the  option  exists  to  have  a  multi-energetic  beam  in  which  case  the  flux  units  must  be  l/(cm 
sec  MeV).  See  electron  energy  flux  spectrum  below. 

2.  Time  resolution  is  the  time  the  charge  deposited  within  the  dielectric  is  allowed  to  move  due 
to  the  electric  fields  before  the  electric  fields  are  recalculated.  It  seems  a  time  step  of  0. 1 
seconds  works  well. 

3.  Spatial  resolution  (cm)  is  the  distance  between  the  nodes  or  bins  that  the  material  is  divided 
into.  Charge  is  considered  to  be  stored  in  each  bin.  The  value  must  be  0.001  cm  or  smaller, 
particularly  for  lower  energies. 

Electron  Energy  Flux  Spectrum 

A  comma-delimited  .csv  file  must  be  in  the  same  folder  as  AFNUMIT2.m  to  specify  the  electron 
energy  flux  spectrum.  The  first  cell  of  the  file  always  is  a  zero.  The  rest  of  the  first  line  contains 
the  energy  (MeV)  of  each  channel  of  flux.  The  following  lines  give  specific  information  about 
the  time  that  a  particular  electron  flux  will  be  incident  and  the  actual  values  of  that  flux.  The  first 
cell  contains  the  time  that  that  particular  electron  energy  flux  spectrum  will  run.  The  remaining 
cells  contain  the  electron  flux  l/(cm"  sr  sec  MeV)  for  each  corresponding  energy  channel.  (An 
example  of  this  file  structure  is  given  in  Table  1  in  the  main  text.)  After  each  row  is  simulated,  a 
line  will  be  plotted  in  all  of  the  output  graphs.  If  more  plot  lines  are  desired  in  order  to  observe 
the  progression  of  the  charging  through  time  (for  the  same  electron  energy  flux  spectrum), 
simply  divide  the  simulation  time  into  the  desired  number  of  lines  and  repeat  the  row  that  many 
times. 

For  mono-energetic  beams,  the  comma-delimited  file  only  contains  two  columns.  The  top  left  is, 
as  always,  zero.  The  top  right  cell  is  the  energy  of  the  beam.  The  first  cell  in  the  second  row  is 
the  time  of  simulation,  as  before,  and  the  electron  flux  goes  in  the  second  column.  Note  that  the 
electron  flux  now  must  have  units  of  l/(cm“  sec).  For  multi-energetic  beams  the  same  format  is 
used  as  for  isotropic,  but  the  electron  flux  must  have  units  of  l/(cm  sec  MeV). 

After  AFNUMT2.m  is  started,  it  will  stop  and  ask  the  user  to  choose  the  electron  energy  flux 
spectrum  .csv  file  from  the  directory. 
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Outputs 


AFNUMIT2.m  will  output  four  plots  in  one  .fig  file: 

1 .  Electric  Field 

2.  Charge  Distribution 

3.  Total  Current 

4.  Dose  Rate 

A  Matlab  .fig  file  can  be  modified  and  reformatted  extensively  within  Matlab,  however,  doing  so 
requires  some  skill  with  Matlab.  Therefore,  AFNUMIT2.m  automatically  creates  four  .csv  files 
corresponding  to  the  four  plots  and  saves  them  in  the  current  folder.  In  addition  to  providing 
precise  values  for  analysis,  comma-delimited  files  can  be  imported  into  most  any  software  and 
modified  as  desired.  Before  editing,  it  is  strongly  recommended  that  these  files  be  moved  to 
another  folder  due  to  the  second  purpose  below. 

The  second  purpose  in  creating  the  .csv  results  files  is  that  they  are  used  by  AFNUMI2.m  when 
it  is  decided  to  later  continue  the  simulation  run.  To  do  so,  the  user  must  choose  [more  =  '  y '  ] 
as  discussed  above.  Note  the  files  must  be  in  the  same  folder  as  AFNUMIT2.m. 
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APPENDIX  B  -  Actual  Input  Code 


The  first  100  lines  of  code  are  reproduced  below  so  that  the  reader  can  easily  identify  the  places 
to  which  the  above  instructions  in  Appendix  A  refer.  These  lines  are  mostly  comments  that 
primarily  explain  how  to  use  AF-NUMIT2,  concisely  repeating  the  above  instructions.  All  user 
inputs  are  included  in  these  lines;  the  only  reason  to  modify  any  other  line  of  code  in 
AFNUMIT2.m  would  be  to  actually  change  the  program. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  AF-NUMIT2 

%  The  AF-NUMIT2  program  calculates  the  one-dimensional  motion  of  electrons 
%  within  a  dielectric  material  that  is  placed  between  two  conducting 
%  surfaces.  The  NUMIT  method  divides  the  material  into  spatial  nodes.  For 
%  each  spatial  node: 

%  1)  The  electric  fields  within  the  dielectric  are  calculated. 

%  2)  The  resulting  movement  of  charge  between  the  nodes  is  calculated. 

%  3)  The  new  distribution  of  charge  at  each  node  is  determined. 

%  The  numerical  calculations  are  iterated  for  each  step  in  time.  Hence  the 
%  name  NUMIT  stands  for  "NUMerical  ITeration."  The  approach  was  developed 
%  by  A.R.  Frederickson  at  the  Air  Force  Cambridge  Research  Labs  in  1974 
%  (AFCRL-TR-74-0582 ) .  For  modeling  of  dielectric  charging  by  electron 
%  bombardment,  NUMIT  requires  the  input  from  algorithms  that  predict  the 
%  deposition  profiles  of  the  incident  electrons  and  their  energy. 

% 

%  AF-NUMIT2  is  a  completely  new  version  of  Frederickson ' s  NUMIT  written  by 
%  B.P.  Beecken  under  the  auspices  of  AFRL/RVBX  and  completed  in  2013.  It 
%  features  re-derived  transport  equations  and  the  ability  to  model  an 
%  isotropic  incident  flux  of  electrons.  The  following  functions  are  part  of 
%  AF-NUMIT2 : 


BeeckenTabFred 

getparameterB 

backscatter 
num  der 


algorithm  predicting  electron  and  energy 
deposition  as  a  function  of  angle 
algorithm  interpolates  a  published  table  to  find 
the  value  of  parameter  B  for  electron  deposition 
a  guess  at  backscatter  of  incident  electrons 
takes  numerical  derivative 


% 

% 

% 

% 

% 

% 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clear  all; 
%close  all; 
%clc; 


%  ensure  no  previously  defined  variables  interfere 
%  close  down  all  open  graphing  windows 
%  clears  main  screen  of  any  previous  input 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%  INPUTS  REQUIRED  FROM  USER  %%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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%%%%%%%%%%%%%%  Return  to  previous  simulation  and  continue? 

%  By  default,  every  simulation  is  new  so  "more"  =  ' n',  which  stands  for  no 
%  or  new. 

%  If  the  user  wishes  to  continue  a  previous  run,  set  "more"  =  ' y '  ,  for  yes. 
%  During  a  continuation,  the  energy  flux  spectrum  can  change,  but  all 
%  material  characteristics  must  be  unchanged, 
more  =  ' n ' ; 

Set  material  characteristics  of  dielectric 


L  =  0.16; 

Z  =  6.3597; 

AW  =  12.3922; 
RHOjyLAT  =  1.42; 
V  =  0; 

gDark  =  7e-22; 
gCoef  =  4.5e-16; 


%  total  thickness  of  dielectric  (cm) 

%  effective  atomic  number  of  dielectric 
%  effective  atomic  weight  of  dielectric 
%  dielectric  density  (g/cmA3) 

%  voltage  difference  between  electrodes  (V) 

%  dark  conductivity  (l/ohm*cm) 
coef  of  radiation  induced  conductivity  ( sec/ohm* cm*rad) 


epsilon  =  3 . 4*8 . 854e-12;  %  permittivity  (epsilon)  of  dielectric  (F/m) 


S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

ooooooooooooooo 


Set  Model 


Parameters  Controling  Resolution 


S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

ooooooooooooooo 


%%%  Set  Isotropic  Incidence  Resolution  of  Model  %%% 

%  Set  number  of  angular  bins:  use  1  for  beams. 

%  For  multi-energetic  beams,  the  flux  must  be  1/ (cmA2*sec*MeV) . 
%  For  mono-energetic  beam,  the  flux  must  be  1/ (cmA2*sec) . 
numberOf AngularBins  =  1; 


%%%  Set  Time  Resolution  of  Model  %%% 

delta_t  =0.1;  %  length  of  each  time  step  (seconds) 

%%%  Set  Spatial  Resolution  of  Model  %%% 

delta_x  =  0.001;  %  distance  between  spatial  nodes  (cm) 


logarithm  =  ' y' 


Energy  Spectrum  &  Incident  Electron 
Input  with  a  CSV  "Data"  file 
Energy  channel  width  linear  or 


Fluxes 


S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

ooooooooooooooo 


S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

oooooooooooooooooooooo 


logarithmic  ('n'  or  ' y ' ) 


O, 

o 

%  INPUT  CSV  FILE  LOCATION:  must  be  in  same  folder  as  AF-NUMIT2 
%  INPUT  CSV  FILE  DESCRIPTION:  The  cell  that  is  the  first  col  and  row  (i.e., 

%  "top  left")  is  ignored,  but  put  a  zero  in  it  because  comma-delimited 
%  files  need  a  place  holder.  The  first  row  of  the  matrix  indicates  the 
%  energy  (MeV)  of  each  channel.  The  second  row  and  all  following  rows 
%  contain  the  electron  flux  in  1/ (cmA2*sr*sec*MeV)  at  the  particular  energy 
%  designated  in  the  first  row,  but  the  first  cell  is  devoted  to  the  time  in 
%  seconds  to  be  simulated  at  that  energy  spectrum. 

%  MULTI-ENERGETIC  BEAMS:  in  this  case  the  input  CSV  file  is  the  same  except 
%  that  the  electron  flux  in  each  energy  channel  must  be  in  1/ (cmA2*sec*MeV) . 
%  No  Steradians !  It  is  a  beam. 

%  MONO-ENERGETIC  BEAMS:  in  this  case  the  input  CSV  file  has  only  2  columns! 

%  In  the  first  row,  as  before,  the  first  cell  is  zero  and  the  second  cell 
%  is  the  energy.  After  that,  the  first  column  is  the  time  (sec)  and  the 
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%  second  column  is  the  electron  flux  in  1/ (cmA2*sec) .  No  steradians.  No 
%  energy  beam  width. 

%  FORMATTING  OF  PLOTS:  after  each  row  (which  corresponds  with  an 
%  energy  spectrum),  NUMIT2  updates  the  plots.  If  more  updates  are  desired 
%  (lines  in  the  plots)  for  the  same  energy  spectrum,  simply  create 
%  identical  energy  spectrum  rows  with  appropriate  times  for  plot  intervals. 
%  OUTPUT  CSV  FILE  LOCATION:  placed  in  same  folder  as  AF-NUMIT2 

S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  TT  Q  T?  P  TTVTPTTTQ  fTiMPT  Th’TTh’Pi  9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 9- 
oooooooooooooooooooooooooo  UOIjI\  -L  IN  Jr  U  -L  O  J_iIL  ±  ILU  oooooooooooooooooooooooooo 
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Directly  Determined  by  User 
Inputs 


APPENDIX  C  -  Important  Variables  in  NUMIT2 


Variable  Name 

Variable 

Type 

Meaning/Purpose 

Source 

L 

number 

Total  thickness  of  dielectric  (cm). 

User  specified. 

Z 

number 

Effective  atomic  number  of 

dielectric. 

User  specified. 

See  instructions  for 
calculating. 

AW 

number 

Effective  atomic  weight  of 
dielectric. 

User  specified. 

See  instructions  for 
calculating. 

gDark 

number 

Dark  conductivity  (g/ohm  cm) 

User  specified. 

gCoef 

number 

Coefficient  of  radiation  induced 
conductivity  (sec/ohm  cm  rad) 

User  specified. 

RHO_MAT 

number 

Dielectric  density  (g/cm3) 

User  specified. 

epsilon 

number 

Permittivity  of  dielectric  (F/m) 

User  specified. 

numberOfAngularBins 

Number 

Resolution  of  isotropic  model. 

User  specified. 

Setting  to  1  models  a  beam. 
Caution:  for  beams  make  sure 

electron  flux  entered  in  Data 
array  is  in  units  of  l/(cm2  sec). 

delta_t 

number 

Length  of  each  time  step  (sec). 

User  specified. 

Tradeoff  accuracy  vs.  run 
time.  If  oscillation  results, 
decrease. 

delta_x 

number 

Distance  (cm)  between  the  spatial 
nodes  at  which  all  values  are 

calculated. 

User  specified. 

Tradeoff  between  spatial 
resolution  and  run  time. 

Data 

array 

User  specified  energy  spectra 
(MeV)  and  corresponding  electron 
fluxes  l/(cm2  sr  sec  MeV).  First  cell 
blank.  First  row  is  the  energy 
channels.  Each  subsequent  row  is 
another  flux.  The  first  column  of 

each  row  is  the  number  of  seconds 
for  that  spectrum  to  be  incident. 
[Note  for  beams  electron  flux  must 
be  in  l/(cm2  sec).] 

User  input  CSV  file. 

spectrumNumber 

number 

Increment  used  to  keep  track 
within  largest  loop  the  energy 
spectrum  being  used. 

Increments  between  1  and  the 
numberOf  EnergySpectra  at 
beginning  of  main  loop. 

numberOf  Energy  Bins 

number 

The  number  of  energy  channels  in 
the  spectrum  (each  with  a  specified 
flux). 

Calculated  from  the  number 
of  columns  in  Data,  less  one. 

TO 

row 

vector 

Energy  of  electrons  (MeV)  for  each 
channel  in  spectrum. 

Each  row  of  Data,  after  the 
first  column. 

numberOf  EnergySpectra 

number 

Number  of  potentially  different 
energy  spectrums  and  the  number 
of  lines  on  final  plot. 

Calculated  from  the  number 
of  rows  in  Data,  less  one. 
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timeForEachEnergySpectrum 

col 

vector 

Time  for  flux  from  a  particular 
energy  spectrum  to  be  incident  on 
dielectric.  It  also  is  time  between 
output  plot  lines. 

First  column  of  Data. 

JJncident 

array 

Incident  electron  flux  [l/(cm2  sr  sec 
MeV)]  for  an  energy  channel.  Each 
row  could  be  a  different  flux,  each 
column  a  different  energy. 

All  but  first  column  and  first 

row  of  Data. 

JJncid 

number 

Incident  electron  current  (A/cm2) 
for  a  specific  energy  bin  within  a 
specific  spectrum. 

Calculated  by  taking  a  cell 
within  the  Data  matrix  and 
multiplying  by  appropriate 
conversion  factors. 

x_value 

col 

vector 

Distance  (cm)  of  each  spatial  node 
from  surface. 

Calculated  by  adding  one 
delta x  for  each  spatial  node. 

x num 

number 

Number  of  spatial  nodes 

Size  of  vector  x value. 

EDEPOS 

row 

vector 

Energy  deposited  (MeV*cm2/g)  per 
electron  in  each  spatial  node. 

Calculated  by  function 
TabFredBeecken  for  each 
energy  and  angle  of  incidence. 

Jofrac 

row 

vector 

Fraction  of  incident  electrons 
passing  through  at  each  node 

Calculated  by  function 
TabFredBeecken  for  each 
energy  and  angle  of  incidence. 

EDEPOS_angles 

array 

Energy  deposited  per  electron  on 
each  node  (row)  for  a  specific  angle 
(col). 

Calculated  from  EDEPOS  by 
taking  the  fraction  of  flux  in 
that  particular  angular  bin. 

Jofrac_angles 

array 

Fraction  of  incident  flux  on  each 
node  (row)  for  a  specific  angle  (col). 

Calculated  from  Jofrac  by 
taking  the  fraction  of  flux  in 
that  particular  angular  bin. 

EDEPOSJso 

col 

vector 

Total  energy  deposited  per  electron 
after  summing  all  angles 

Calculated  by  summing 
EDEPOS_angles.  Could  be 
replaced  at  this  point  by 
different  algorithm  source. 

Jofracjso 

col 

vector 

Total  fraction  of  incident  electrons 
after  summing  all  angles 

Calculated  by  summing 
Jofrac_angles.  Could  be 
replaced  at  this  point  by 
different  algorithm  source. 

rho_iso 

col 

vector 

Initial  fraction  of  incident  electrons 
deposited  per  unit  time  (C/(cm3 
sec)) 

Derivative  of  Jofracjso.  Could 
be  replaced  at  this  point  by 
different  algorithm  source. 

Dose_SumAngles 

col 

vector 

Dose  (rad/sec)  at  each  spatial  node. 

Calculated  from  EDEPOSJso 
for  each  JJncid. 

Jo_SumAngles 

col 

vector 

Initial  current  (A/cm2)  at  each 
spatial  node. 

Calculated  from  Jofracjso  for 
each  JJncid. 

rho_SumAngles 

col 

vector 

Initial  charge  deposited  (C/cm3)  at 
each  spatial  node. 

Calculated  from  rhojso  for 
each  JJncid  and  multiplying 
by  time  increment.  (Continuity 
Eq.) 

Dose_energies 

array 

Net  initial  dose  rate  for  each  energy 
channel  (col)  and  at  each  spatial 

Found  by  summing  along  the 
rows  of  Dose_angles  for  each 
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Outputs  of  Model 


node  (row). 

energy  channel. 

Jo_energies 

array 

Net  initial  current  (A/cm2)  for  each 
energy  channels  (col)  and  at  each 
spatial  node  (row). 

Found  by  summing  along  the 
rows  of  Jo_angles  for  each 
energy  channel. 

rho_energies 

array 

Net  charge  density  (C/cm3)  for  each 
energy  channel  (col)  and  at  each 
spatial  node  (row). 

Found  by  summing  along  the 
rows  of  rho_angles  for  each 
energy  channel. 

Dose 

col 

vector 

Total  dose  rate  for  all  angles  and 
energies  in  one  spectrum. 

Sum  across  rows  of 
Dose_energies. 

Jo 

col 

vector 

Total  initial  current  for  all  angles 
and  energies  in  one  spectrum. 

Sum  across  rows  of 

Jo_energies. 

rho 

array 
(2  cols) 

Net  charge  density  (C/cm3)  at  each 
node,  including  previous  charge 
depositions.  Each  row  is  node,  and 
each  column  is  a  time  step. 

Sum  across  rows  of 
rho_energies  and  add  to 
previous  depositions  and/or 
results  of  charge  transport 
(i.e.,  the  time  loop). 

E_Total 

col 

vector 

Total  Electric  field  at  each  spatial 
node  within  dielectric. 

Calculated  in  time  loop  from 
charge  within  dielectric  and 
induced  on  electrodes. 

J 

col 

vector 

Total  current  at  each  spatial  node 
within  dielectric. 

Calculated  in  time  loop  using 
Fowler  model. 

E_f  i  n  a  1 

array 

Electric  field  (V/cm)  after  each 
energy  spectrum  (col)  to  be  plotted. 
Each  row  is  a  spatial  node. 

A  col  is  appended  for  each 
E_Total  calculated  for  each 
spectrum. 

CD 

T3 

O 

2 

s— 

J_f  i  n  a  1 

array 

Current  (A/cm2)  after  each  energy 
spectrum  (col)  to  be  plotted.  Each 
row  is  a  spatial  node. 

A  col  is  appended  for  each  J 
calculated  for  each  spectrum. 

O 

to 

+■» 

3 

Q. 

+-» 

o 

rho_final 

array 

Cumulative  charge  density  (C/cm3) 
after  each  energy  spectrum  (col)  to 
be  plotted.  Each  row  is  a  spatial 
node. 

A  col  is  added  for  each  rho 
calculated  for  each  spectrum. 

Dose_final 

array 

Net  dose  rate  (rad/sec)  after  each 
energy  spectrum  (col)  to  be  plotted. 
Each  row  is  a  spatial  node. 

A  col  is  appended  for  each 

Dose  calculated  for  each 

spectrum. 
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APPENDIX  D  -  Calculating  Effective  Atomic  Number  and  Weight  for  Compounds 

Tabata  defines  the  effective  atomic  number  as 


and  the  effective  atomic  weight  as 


^eff  — 


•^eff 


St/ 


lA; 


Tabata  says,  “f  is  fraction  by  weight  of  the  constituent  elements.”  Following  Tabata’s  reference 
to  Rao,  we  find 


r  _  rijAij 

'l  M 

Here  n,  is  the  number  of  atoms,  A '  is  the  atomic  weight,  and  M  is  the  molecular  weight  of  the 
compound  is  given  by 


M  =  XjrijA'j  . 

For  example,  the  compound  Kapton  C22O5N2H10  has  molecular  weight  M  calculated  by 
M  =  22(12.011)  +  5(15.9994)  +  2(14.0067)  +  10(1.00794)  =  382.332 


The  fraction  by  weight  of  the  constituent  elements  are  then: 


fc 


fo 


22(12.011) 

=  0.6911 

382.332 

/tv  — 

5(15.9994) 

=  0.2092 

382.332 

fH  = 

2(14.0067) 

382.332 


0.0733 


10(1.00794) 

382.332 


0.0264 


Finally  then 


Zeft  =  fc(  6)  +  M  8)  +  fN(  7)  +  fH( 1)  =  6.3597 
and 


6.3597 

fci  6/12)  +  /0(8/12)  +  /w(  7/14)  +  fH(  1/1) 


12.3922 
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